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CHAPTER I 


INTRODUCTION 

1.1 MOTIVATION AND BACKGROUND. 

The computer simulation of electromagnetic (EM) field problems 
has received considerable attention during the past thirty years. This is a 
direct result of the realization among scientists and engineers that the com- 
puter could be exploited to solve EM problems that had previously been too 
complex or cumbersome to be treated by established analytical techmques. 
Before the advent of computer technology, solutions of EM field problems 
■were often limited to relatively simple problems or geometries. Generally, 
the intractability of standard analytical techniques and the nonlinearity of 
solutions of EM problems is largely a result of the coupling between electric 
and magnetic fields. With the aid of computers, the inherent nonlinearity 
of EM field problems became more manageable and more complex EM sys- 
tems problems could be and were solved. When problem solving capabilities 
advanced, so did computer and electrical engineering technology. These new 
technologies were then increasingly used in aerospace applications. These 
applications were in control and guidance as well as more efiBcient servos, mo- 
tors, generators and electronic sensing and surveillance gear. More recently, 
high-temperature superconductors (HTS) have been discovered. These com- 
posite materials are presently the subject of intensive experimental research 
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and are expected to have a major impact in space propiilsion, digital com- 
puting, power systems and communications in the next century. 

The present work is part of a research program for the numerical 
simulation of EM/mechanical systems that involve superconductors. The 
point of departure from previous works is the use of finite elements based 
upon a four-potential variational principle to predict desired EM quantities. 

The simulation involves the interaction of the following four compo- 
nents: 

(1) Thermal Fields: temperature and heat fluxes. 

(2) Electromagnetic Fields: electric and magnetic field strengths and fiuxes, 
currents and charges. 

(3) Quantum Mechanics: the constitutive behavior of the superconducting 
system is governed by queintum mechanical effects. Particularly impor- 
tant is the superconducting phase change, governed by phenomena at 
the quantum level, and triggered by thermal, mechcinical and EM field 
energy levels. 

All three components can be treated by the finite element method. 
This treatment produces the spatial discretization of the continuum into 
mechanical, thermal, quantum mechanical and electromagnetic meshes of a 
finite number of degrees of freedom. The fimte element discretization may 
be developed in two ways: 

(1) Simultaneous Treatment. The whole problem is treated as an indivisi- 
ble whole. The four meshes noted above become tightly coupled, with 
common nodes and elements. 
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(2) Staged Treatment. The mechajaical, thermal and electromagnetic com- 
ponents of the problem are treated separately. Finite element meshes 
for these components may be developed separately. Coupling effects 
are viewed as information that has to be transferred between these four 
meshes. 

The present research follows the staged treatment. More specifically, 
we develop finite element models for the fields in isolation, and then treat 
coupling effects as interaction forces between these models. This “divide 
and conquer” strategy is ingrained in the partitioned treatment of coupled 
problems [2,3], which offers significant advantages in terms of computational 
efficiency and software modularity. Another advantage relates to the way 
research into complex problems can be made more productive. It centers 
on the observation that some aspects of the problem are either better un- 
derstood or less physically relevant than others. These aspects may then be 
temporarily left alone while efforts are concentrated on the less developed 
and/or more physically important aspects. The staged treatment is better 
suited to this approach. Of the foxir components listed previously, the last 
two are less developed in a modeling and computational sense. 

Mechanical elements for this research have been derived using general 
variational principles that decouple the element boimdary from the interior, 
thus providing efficient ways to work out coupling with non-mechamcal fields. 
The point of departure was the previous research into the free-formulation 
vzfriational principles presented by Felippa [4]. A more general formulation 
for the mechanical elements, which includes the assumed naturcd deviatoric 
strain formulation was established and reported in Refs. [5, 6, 7,8]. New repre- 
sentations of thermal fields have not been addressed zis standard formulations 
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axe considered adequate for the coupled-field phases of this research. How- 
ever, research in thermomechanical interactions supported by this program 
has resulted in the construction of robust and efiicient staggered solution 
procedures [9]. 

The development of EM finite elements to date has not received the 
same degree of attention given to mechanical and thermal elements. Part of 
the re£ison is the widespread use of analytical and semianalytical methods 
in electricad engineering. These methods have been highly refined for spe- 
cialized but important problems such as circuits and wave guides. Thus the 
advantages of finite elements in terms of generality have not been enough to 
coxmterweight established techniques. Much of the EM fimte element work to 
date has been done in England and is well described in the surveys by Davies 
[10] and Trowbridge [11]. The general impression conveyed by these surveys 
is one of an unsettled subject, reminiscent of the early period (1960-1970) 
of finite elements in structural mechamcs. A great number of formulations 
that combine fiux, intensity, emd scalar potentials are described with the 
recommended choice varying according to the application, medium involved 
(polarizable, dielectric, semiconductors, etc.), number of spatial dimensions, 
time-dependent characteristics (static, quasi-static, harmonic, or transient), 
as well as other factors of lesser importance. The possibility of a general 
variationzil formiilation has not been recogmzed. 
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1.2 T^EVTEW OF EXISTING TECHNIQUES 

As mentioned previously, the computer simulation and modeling of 
EM field problems is presently an unsettled subject especially in nonlinear 
problems. A rich variety of mathematical techniques have been used to solve 
these complex problems. Some of these techniques involve using integral 
transforms to find a solution, while other techniques yield solutions to the 
integral or differential EM field equations that contain Bessel, Airy, Gamma, 
and Legendre functions [12]. A common method of computer implementa- 
tion involves taking the analytical representation of the solution to a problem 
and making a numerical approximation to that solution. In time-independent 
problems the implementation may take the form of discretizing the analytical 
differential or integral EM field equations over the system’s spatial dimen- 
sions [13]. Linear time-dependent problems may be transformed to Fourier 
or Laplace space, solved, and then converted back to the real time domain. 
The computer is simply used to make good approximations to an integral 
which is an analytical solution to the problem, but for which no closed form 
solution of the integral exists. While these methods are effective for specific 
problems, they axe rarely of a general enough nature that they can be used 
on most EM system problems. A recogmtion of the interest in and the need 
for more generalized computer solution techniques for EM field problems led 
to the first COMPUMAG series of conferences in 1976 [11, p. 506]. 

Prior to this time, few finite element techniques existed but the power 
of more generalized schemes were demonstrated in finite difference codings 
that used the differential forms of Maxwell’s field equations. Usually the 
conventional field quantities were replaced by potentials and the resultant 
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EM field equations were discretized over space [14, pp. 101-105;15,16]. Fi- 
nite difference schemes generally axe not as amenable to Neumann boundary 
conditions or an easy change to a higher variational order as finite element 
methods. Because of the prevelance of Neuma nn bovmdary conditions in EM 
field problems, especially when conventional field quantities are replaced by 
a potential formulation, as well as difficulties associated with a change in the 
variational order of variables when finite difference methods axe used, finite 
difference techniques axe reirely used for the spatial discretization of EM field 
equations. 

Maxwell’s EM field equations may be recast in a potential formula- 
tion. This reduces the number of independent variables for the electric field 
E from three to one through the substitution E = -V$, where $ is the 
electrostatic potential. The reformulation of the magnetic field is more com- 
plicated. In free space, the magnetic field B can be defined as the negative 
gradient of the magnetostatic potential <p (t.c., B= — This substitu- 
tion reduces the number of independent variables from three to one for the 
magnetic field but this potential is neither single valued nor defined in a 
conductor that is carrying a steady current [14, p.l39]. Another reformu- 
lation substitutes the curl of the magnetic vector potential A for B (i.e., 
V X A = B). Although this formulation does not reduce the number of inde- 
pendent variables in an EM field problem, it does require that the solution 
of A be C° continuous across material interfaces, thus simplifying finite ele- 
ment development. Formulations that use the B field as a primary variable 
are not required to be C° continuous across material boundaries. In spite 
of the difficulties presented by a discontinuous variable, the majority of EM 
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field finite element formulations to date are based on the original EM fields, 
e.g., see Refs. [17,18]. 

Some researchers have also experimented with magnetic vector po- 
tential bzised finite elements [11]. These formulations use a Galerkin weighted 
residual method applied to the strong form of the EM field equations. The 
drawback with this approach is that the xiniqueness of a numerical solution 
is questionable because the divergence of A is not specified. A variational 
approach based upon A can easily overcome this difficulty by specifying a 
function or gauge for the divergence of A, weighting it by a Lagrangian mul- 
tiplier, and augmenting it to the energy functional of the EM system. The 
only requirement on the choice of gauge is that the Euler equations of the 
weighted gauge choice equal zero. Another statement of this requirement is 
that the augmented energy functional should differ from the EM field en- 
ergy functional by a constaint [19, p. 36]. In fact, by an appropriate use 
of the Loreniz gauge, the Lagrangian, or energy functional, of the EM field 
equations can be used to perform a canonical transformation to produce the 
Hamiltonian of the system [20, pp. 72-91]. The only EM finite elements 
that use the approach of energy functionals augmented by a weighted gauge 
equation are the ones presented in this work. 

As mentioned on the previous page, the magnetic scalar potential 
can be used to calculate the B field in free space and reduce the number of 
independent variables from three to one. To increase computational speed 
and reduce memory allocation, Trowbridge [11] has coupled A with to 
produce a new independent variable vector quantity R. R requires that three 
variables be solved in a conducting media zind only one variable be solved in 
free space. This method has drawbacks, specifically that R and <p are not 
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unique, and that in the interior of conductors cancellation problems arise 
that can give erroneous values for B [11, pp. 521-525]. 

To model EM fields in a superconductor, another field variable must 
be included: the wave order par 2 imeter rj}. This function CEin be complex and 
the absolute value of ip times its complex conjugate (\rpip*\ = |^P) is defined 
as the number density of superconducting electron pairs. This new variable 
accounts for the quantum mechamcal effects that appear in the interior of a 
superconductor. These quantum effects change the value of the B field and 
ciirrent density vector j within a superconduct^' “ 

A widely used mathematical model that describes quantum and EM 
interactions within a superconductor are the Ginzburg- Landau equations [21, 
p. 104] . These equations reduce to Maxwell’s equations, the same equa- 
tions that govern EM fields in normal conductors and in vacuum. The 
Ginzbmg-Landau equations are derived using variational principles and re- 
quire a unique gauge choice to ensure that a superconducting current can 
only exist in a conductor as physics demands. The gauge choice used in the 
present work is called the London gauge and is equivalent to the Lorentz 
gauge for magnetostatic problems. The Ginzburg-Landau equations also 
contain A explicitly as well as the vector curl of A. To model superconduc- 
tors numerically, the optimal choice of independent variables is to use A. 
The use of a field based formulation requires the numerical integration of 
B to remove terms in the Ginzburg-Landau equations that contain A. This 
integration can easily become the source of additional numerical error, an 
error that is not present when the choice of independent variable is A. 

A finite difference formulation of the Ginzburg-Landau equations has 
been developed that producs reasonable results [15,16]. The formulation uses 
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A as a primary variable, but thermal elFects eure neglected when this model is 
in the normal state. This formulation also suffers from the previously men- 
tioned drawbacks of finite difference methods in the treatment of arbitrary 
geometries. 

1.3 CONTENT. 

The objective of this thesis is to develop EM finite elements for type I 
anH II superconductors based upon a gauged four-potential variational prin- 
ciple. At present, the physics of high temperature superconductors (HTS) are 
not well enough imderstood to permit the development of zm adequate math- 
ematical model. The last elements developed in this work include thermal 
coupling, but are magnetostatic. This restriction is motivated by the fact 
that the time-independent problem exhibits strong nonlinearities; further- 
more, no completely satisfactory mathematical model has been developed 
for the time-dependent case [21, p. 273]. The highly nonlinear nature of 
the problem is the result of a boundary layer effect exhibited at a supercon- 
ductor/normal conductor or superconductor /vacuum interface. Extremely 
strong gradients of the independent variables V’l -A-j B, and j are present in 
this regime. These gradients bring about serious numerical diflB.culties, the 
most important ones being a highly ill-conditioned system of incremental 
equations and the need for specialized mesh discretization. The final super- 
conducting finite element developed is of a general enough nature that it 
works equally well in both the boundary layer and the bulk of the super- 
conductor. Unlike the previously mentioned field based formulations, this 
element requires no special treatment for material interfaces, in particular, 
the superconductor /vacuum interface. 


iiiiiL 
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The derivation of edl of the EM finite elements in this thesis are based 
upon a four-potential variational formulation that uses the four-potential as 
the primary variable. The electric field is represented by a scalar potential 
and the magnetic field by a vector potential. When the superconductor is 
modeled, the electric field sceJar potential is dropped, because it does not 
couple with the magnetic field in the magnetostatic case. The modulus and 
phase of ip are then added as new independent variables. The formulation of 
the four-potential variational principle proceeds along lines previously devel- 
oped for the acoustic fluid problem |22,23]. The appropriate gauge normal- 
ization is incorporated in the variationzd (weak) form through the adjimction 
of a Lagrange multiplier field. 

The main advantages of developing finite elements using a potential 
based variational formulation in contrast to using existing EM numerical 
techniques are summarized as follows. 

(1) Interface discontinuities are automatically taken care of without any 
special intervention. 

(2) No approximations are invoked a priori since the general Maxwell equa- 
tions are used. 

(3) The number of degrees of freedom per finite element node is kept modest 
as the problem dimensionality increases. 

(4) Higher order and hybrid elements are more eeisily aocomodated. 

(5) The Ginzburg-Landau equations naturally possess A as an independent 
variable; possibilities for errors from an additional numerical integration 
are removed. 

(6) A generalized formulation that posesses a broad range of applicability. 
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REMARK 1.3.1 

An interesting byproduct of this formulation is that with minor modifications, it 
can be used to describe the physics of a superfluid. See Ref. [20], pp. 152-158. 

1.3.1 FINITE ELEMENTS. 

A total of eight finite elements were developed in the course of the 
the author’s research. Seven of these are based upon the four-potential varia- 
tional principle, and the last is a thermal conduction element developed from 
a different variational principle according to Ref. [24]. They are in order of 
development: 

(1) a one- dimensional Coupled Linear Electric and Magnetic field (CLEMlD) 
finite element 

(2) a two-dimensional axisymmetric Coupled Linear Electric and Magnetic 
field (CLEM2D) finite element 

(3) a one-dimensional Coupled Linear Electric and Magnetic field INFinite 
(CLEMINF) finite element 

(4) a one- dimensional CUrrent Predicting Linear Electromagnetic (CU- 
PLEID) finite element 

(5) a one- dimensional Superconducting Thermal, Electromagnetic and Phase 
coupled (STEPID) finite element 

(6) a one-dimensional Superconducting ThErmAl, and eLectromagnetic 
field (STEALID) finite element 

(7) a one-dimensional LINear Thermal conduction (LINTlD) finite element 

(8) a one-dimensional Linear Electromagnetic and Thermally coupled (LETlD) 
finite element 
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Elements (1), (2) eind (3) predict only electric and magnetic fields. Element 
(3) was developed as a term project, but has limited practical usage except 
for the development of an EM finite element that is time-dependent. Ele- 
ments (4), (5), (6), and (8) can predict EM fields, but also have the ability 
to predict the current density distribution, j, given the scalar input I, the 
total ciirrent. Element (6) is not presented here, because it can easily be 
derived from element (5) by constraining the variable \x})\ to be a constant. 
This formulation is known as the the London formulation for superconduc- 
tors. This element was developed solely for the purpose of troubleshooting 
element(5) [25]. Element (5) also predicts the quantum mechanical quantity 
|^|. It also contains two thermally dependent material parameters. These 
two parameters couple the superconductor to thermal fields. Element (7) 
was constructed to predict the temperature distribution within the conduc- 
tor. Element(S) can predict j and EM fields, but is coupled to thermal fields 

by the electrical resistivity, a;. 

REMARK 1 . 3.2 

Appropriate changes to the Ginzburg- Landau theory and finite element formulation 
for the construction of element (6) are listed in this thesis. Results for (6) are 
deleted as they are not as accurate as the results obtained from the STEP ID finite 
element which is based upon the complete Ginzburg-Landau theory where V’ is 
allowed to vary. 

1.3.2 mSSERTATION OUTLINE. 

The dissertation is organized as follows. Chapter II is devoted to 
a review of basic EM theory, and the development of four-potential theory. 
Variational functionals for two cases where the current density vector j is 
known are also discussed. Chapter III is devoted to the development of vari- 
ational functionals for conductors where j is imdetermined. In this chapter. 
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functionals for the normal and superconducting states of a conductor are pre- 
sented. Chapter IV introduces the variational functionals necessary for the 
time-independent heat conduction and convection problems. Some general- 
ized solutions for one-dimensional conductors are also presented here. This 
chapter also includes formulas that express the values of a conductor’s EM 
material properties as a function of the temperature 'T . Accurate numerical 
approximations for the values of these material properties are also developed. 

The first four chapters outlined above comprise the first step in the 
development of EM finite elements that can model the quantum and ther- 
mal effects that appear within a superconducting material. The main goal of 
these chapters is to develop variational functionals that are later discretized 
to produce finite elements. These fimte elements are then used to analyze the 
thermal, quanttim and electromagnetic properties of a conductor for some 
specific EM field problems. The following seven chapters are devoted to de- 
veloping finite elements and solving those specific EM field problems. Where 
an analytical solution to the field problem exists, it is presented in that chap- 
ter. ff special numerical procedures are necessary for the solution of the field 
problem, the procedures are also discussed in that chapter. Chapters V and 
VI deal with one and two-dimensional axisymmetric EM field problems re- 
spectively, where the current density vector j is known and the conductor 
remains in the normal state. Chapter VII presents the finite element solu- 
tion of a one-dimensional axisymmetric conductor in its normal state where 
the current density vector j is unknown. Chapter VIII is concerned with 
finding the values of EM fields within a one-dimensional time-independent 
axisymmetric superconductor. Chapter DC develops a one-dimensional heat 
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conduction finite element. This element is employed with a modified ver- 
sion of the element of Chapter VII to solve the coupled problem of a one- 
dimensional axisymmetric conductor that is subjected to a varying thermal 
loeid. Chapter X employs appropriately modified versions of the elements 
of Chapter IX to solve the coupled EM-thermal system where the electric 
current through a one- dimensional axisymmetric wire is varied. Chapter 
XI models the complete quantum, themal and EM field problem for a one- 
dimensional axisymmetric wire. The temperature T and the electric current 
are allowed to vary, but the wire is also allowed to change its quantum state 
and be either a normal conductor or a superconductor. 

The last chapter, Chapter XII, contains a broad summary of the 
dissertation. This chapter highlights some of the more important aspects of 
the variational methods used here. It concludes the dissertation with a small 
section on new research directions that the thesis research has suggested. 


CHAPTER II 


EM AND FOUR-POTENTIAL THEORY 


2.1 ELECTROMAONETIC FIELD EQUATIONS. 


2.1.1 THE MAXWELL EQUATIONS. 

The original Maxwell equations (1873) involve four three- vector 
quantities; B, D, E, and H. Vectors E and H represents the electric and 
magnetic field strengths, respectively, whereas D and B represent the electric 
and magnetic fluxes, respectively. All of these are three-vector quantities, 
that is, vector fields in three-dimensional space (e.g., in Cartesian space, 
xi =x, X 2 = y, xz = z): 

( D^^ ( ( Hi ^ 

( 2 . 1 . 1 ) 

Other quantities are the electric current 3-vector j and the electric charge 
density p (a scalar). 

With this notation, and using superposed dots to denote differenti- 
ation with respect to time t, Maxwell equations can be stated as 


E= < 

E, 

> D= ^ 

[g] 

> B= < 


\ K= < 

H2 / 


UsJ 




[BsJ 

1 

i^3 J 


B -1- VxE = 0 

VxH-D =j 

VD = p 

V-B = 0 


( 2 . 1 . 2 ) 


The first and second equation are also known as Faraday’s and Ampere- 
Maxwell laws, respectively. 
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The system (2.1.2) supplies a total of eight partial differential equa- 
tions, which as stated are independent of the properties of the underlying 

medium. 

REMARK 2.1.1 

Some authors, for example, Eyges [26], include Air factors and the speed of light c 
in the Maxwell equations. Other textbooks, e.g. [27, 28], follow Heaviside’s advice 
in using technical units that eliminate such confusing factors. 


nONSTTTTTTTVE EQUATION 


The field intensities E and H and the corresponding flux densities 
D and B are not independent but are connected by the electromagnetic 
constitutive equations. For ein electromagnetically isotropic, non-polarized 
material the equations are 

B = fxB. D = eE I (2.1.3) 

where fj, and e are the permeability and permitivity, respectively, of the ma- 
terial. These coefficients are functions of position but (for static or harmonic 
fields) do not depend on time. In the general case of a non-isotropic mate- 
rial both [1 and e become tensors. Even in isotropic media /i in general is 
a complicated function of H; in ferromagnetic materials it depends on the 
previous history (hysteresis effect). 

In free space fi = fio and e — eo? which are connected by the relation 



(2.1.4) 


where Cq is the speed of light in a free vacuum. In rationalized MKS units, 


Co ~ 3.10® m/sec and 


/io =?= 47t X 10"^ henry/m, cq = Hq ^ = (367t) ^ x 10 sec2/(henry • m) 

(2.1.5) 
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The condition /z w /xq holds well for most practical purposes in such media 
as air and copper; in fact flair — 1-0000004/zo and Hcopper = .99999^0- 

The electric field strength E is further related to the current density 
j by Ohm’s law: 

j = <tE (2.1.6) 

where a is the conductivity of the material. Again for a non-isotropic mate- 
rial a is genereiUy a tensor which may also contain real and imaginary com- 
ponents; in which case the above relation becomes the generalized Ohm’s 
law. For good conductors <r > > e; for bad conductors a < < e. In free space, 
<7 = 0 . 

2.1.3 MAXWELL EQUATIONS IN TERMS OF E AND B. 

To pass to the four-potential considered in this work it is convenient 
to express Maxwell’s e<iuations in terms of the electric field strength E and 
the magnetic flux B. In fact this is the pair most frequently used in elec- 
tromagnetic work that involve arbitrary media. On eliminating D and H 
through the constitutive equations (2.1.3), we obtain 


B -1- VxE = 0 

VxB — fieE = fij 

V • E = p/e 

< 

II 

o 


The second equation assumes that e is independent of time; otherwise eE = 

edE/dt shovdd be replaced by d{eE)ldt. In charge-free vacuum the equations 

reduce to ^ 

B-l-VxE = 0 VxB--%E = 0 

VE = 0 V-B = 0 


( 2 . 1 . 8 ) 
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2.1.4 THE ELECTROMAGNETIC POTENTIALS. 

The electric scalar potential $ and the magnetic vector potential A 
are introduced by the definitions 


E = -V$ — A B = VxA 


(2.1.9) 


This definition satisfies the two homogeneous Maxwell equations in (2.1.7). 
The definition of A leaves its divergence V • A arbitrary. We shall use the 


Lorentz gauge [29] 

V • A H- fie^ = 0 


( 2 . 1 . 10 ) 


With this choice the two non-homogeneous Maxwell equations written in 
terms of $ and A separate into the wave equations 


- /le# = -p/« A - /xeA = -w (2.1.11) 


2.2 THE ELECTHOMAGNETIC FOUR-POTENTIAL. 

Maxwell’s equations can be presented in a compact manner (a form 
compatible with special relativity) in the four-dimensional spacetime defined 
by the coordinates 

X\ =X, I2 = y, ^3 = ^) ^4 = ict (2.2.1) 

where Xi,X 2 ,xz are spatial Cartesian coordinates, = —1 is the imaginary 
unit, and c = is the speed of E M waves in the medium under con- 

sideration. In the sequel Roman subscripts will consistently go from 1 to 4 
and the summation convention over repeated indices is used unless otherwise 

stated. 
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2.2.1 THE FIELD STRENGTH TENSOR. 


The unification can be expressed most conveniently in terms of the 
field- strength tensor F, which is a four-dimensional antisymmetric tensor 
constructed from the components of E and B as follows; 


■ 0 

Fi2 

Fi3 

Fn' 


■ 0 

cBz 

—cBz 

-iEi ' 

-Fn 

0 

F 23 

F 24 

def 

—cBz 

0 

cB\ 

-iEz 

-Fu 

—F 23 

0 

F 34 

= 7 

CB 2 

— cJ3i 

0 

-iEz 

-Fu 

—F 23 

-F3A 

0 , 


iE\ 

iEz 

iEz 

0 


( 2 . 2 . 2 ) 


Here 7 is an adjustment factor to be determined later. Similarly, we can 
introduce the four-current vector J as 


fJil 


' cm ' 


' m '' 

J2 

J3 

del 

► = 7 < 

Cfljz 

cm 

> = 7c ^ 

m 

m 

U4J 




iy/php J 


(2.2.3) 


Then, for arbitrary 7, the non-hompgeneous Maxwell equations, namely 
VxB — /xeE = and V • E = p/e, may be presented in the compact 
“continuity” form (the covariant form of these two equations): 


^ = Ji (2.2.4) 

dxk 


The other two Maxwell equations, V • B = 0 and VxE -t- B = 0, can be 


presented as 


dFik dFmi dFkm 

dxm dxk dxi 


(2.2.5) 


where the index triplet {i,k,m) takes on the values (1,2,3), (4,2,3), (4,3,1) 


and (4,1,2). 
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2.2.2 THE FOUR-POTENTIAL. 

The EM “four-potential” is a four-vector whose components are 
constructed with the electric and magnetic potential components of A and 




cAi 


, 1 <^2 I del I cA2 

<^ = ^^3 f = ic.43 


( 2 . 2 . 6 ) 


<^4 J 

It may then be verified that F can be expressed as the four-curl of <^, that is 

(2.2.7) 

dxi dXk' ^ 

or in more detail emd using commas to abbreviate partizJ derivatives: 


F = 


0 ^2,1 — 4>3,1 — ^1,3 4>4,1 — 4>\A 

<f>l,2 - 4>2,1 0 4z, 2 — <f>2,Z <i>4,2 ~ <l>2,4 

4>1,3 - <i>3,l <^2,3 - <f>3,2 0 ~ 4>3,4 

L<^1,4 — <f>4,l <i>2,4 — <l>4,2 <i>3,4 ~ <f>4,3 0 


( 2 . 2 . 8 ) 


2.2.3 THE UNGAUG ED LA<^NG!AN. 

With these definitions, the basic Lagrangian of electromagnetism can 
be stated as 


L-^FikFik -47 


(2.2.9) 


= ^7^ (c^B^ — F^) — ^(j\Ai d- j2-^2 + ~ P^) 


in which 


= Bl, = E^E = -|- F| -h F| (2.2.10) 


Comparing the first term with the magnetic and electric energy densities 

[26,27,28] 




Ue = |D^E = ieF^, 


( 2 . 2 . 11 ) 
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we must have — 7^/(^e) = l//i, from which 

7 = Vi (2.2.12) 

Consequently, the required Lagrangian is 

L — ~ ~ d" j2-^2 + j3-^3 ~ /^^)> (2.2.13) 

2fi ‘‘ 

The associated variational form is 

R= r f LdVdt (2.2.14) 

Jto Jv 

where V is the integration volume considered in the analysis. In theory V 
extends over the whole space, but in the numerical simulation the integration 
is truncated at a known boundary or special devices are used to treat the 

decay behavior at infinity. 

REMARK 2.2.1 

Lanczos [30] presents this Lagrangian for free space, but the expression (2.2.13) 
for an arbitrary material was found in none of the textbooks on electromagnetism 
listed in the References. 

2.2.4 THE GAUGED LAGRANGIAN. 

If the fields A and $ to be inserted into L do not satisfy the Lorentz 
gauge relation (2.1.10) a priori, this condition has to be imposed as a con- 
straint using a Lagrange multiplier field Xg{xi), leading to the modified or 
“gauged” Lagrangian: 


Lg = L Aj(V • A -|- 


(2.2.15) 
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2.2.5 THE FQUR-FTFXD EQUATIONS. 

On setting the variation of the functional (2.2.15) to zero we re- 
cover the field equations (2.2.4) and (2.2.5), as well as the gauge constraint 

(2.1.10) as Euler-Lagrange equations. Taking the divergence of both sides of 
(2.2.4) and observing that F is an antisymmetric tensor so that its divergence 
vanishes we get 

|^=cKV-j + ^) = 0 (2.2.16) 

The vanishing term in parenthesis is the equation of continuity, which ex- 
presses the law of conservation of charge. The Lorentz gauge condition 

(2.1.10) may be stated as Finally, the potential wave equations 

(2.1.11) may be expressed in compact form as 

n<i>i = -Ji (2.2.17) 


where □ denotes the “fotir- wave-operator”, also called the D’Alembertian: 


^ d.. s‘ , S‘ , 3^ 

dxkdxk dx\ dx^ dx^ c^df^ 


(2.2.18) 


Hence each component of the four-potential <j> satisfies an inhomogeneous 
wave equation. In free space, = 0 and each component satisfies the homo- 
geneous wave equation. 

The following sections of this chapter are devoted to derivations of 
the appropriate expression for Xj for selected cases. The first variation of R 
with respect to the independent variables is also taken. With few exceptons, 
the solutions of the independent variables $ is not determined. The variation 
is performed primarily to determine the natural boimdary conditions of each 
test case for the eventual extension of the four-potential method to finite 
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element analysis. The variation is also performed to confirm the validity of 
the four-potential method as an analytical tool by directly comparing the re- 
sultant Euler equations minus the Lagrange multiplier terms with Maxwell s 
field equations. 

2.3 THE ONE-DTMENSIONAL AXISYMMETRIC CASE. 

The simplest application for the four-potential variational principle 
is to an infinitely long, straight conductor of circular cross section which car- 
ries a known, time-independent uniform current in the longitudinal direction 
(Fig. 2.1.). To take advantage of the axisymmetric geometry a cylindrical 
coordinate system is chosen with the wire centerline as the longitudinal z- 
axis. The vector components in the cylindrical coordinate directions r , B and 
z are denoted by 

El = Ari Br, Er in the r (radial) direction, 

A 2 , B 2 , E 2 = Ae, Be, Ee in the 6 (circumferential) direction, 

Az, Bz, Ez = A^, Bt, E 2 in the z (longitudinal) direction. 

The first step in solving for the fields is to express the gauged La- 
grangian 

L, = ^B^ - \eE^ - (j^A - p$) -h A,(V • A -h /xe$), (2.3.1) 

in terms of the potentials written in cylindrical coordinates. 
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For we use (2.1.9), (2.2.10) and the cylindrical-coordinate curl formulas 
to get 





/a Ar dA,y fl djrA,) 
dz ) \ dr ) \r 


IdArV 

r do J 


(2.3.2) 


For E'^ we use (2.1.9) and the cylindrical-coordinate gradient formulas to 
produce 


E, 


Er 


E={E2} = {Ee)=-{ 


Ez 


E. 


1 5$ , j ^ 

^ I i 


(2.3.3) 


so that (2.2.10) becomes 


= E-E = (§ + ^)^ (if + ^)' + (^ + ^)^ (2-3.4) 


For the Lorentz gauge we use the cylindrical-coordinate divergence formula 
to get 

„ . -1 d(rAr) 1 dAs dA^ . • 

V • A -f- fj.e^ = 1 4- h (2.3.5) 

r dr r du dz 

The electromagnetic fields, for the one-dimensional czise, only vary in the 
radial (r) direction and any partials with respect to 6 and z vanish. In the 
time-independent case, all partials with respect to t also vanish. With no 
static charge density, /? = 0, and with only a longitudinal current, the single 
non- vanishing component of j is jz. The constitutive relation (2.1.6) can be 
used to remove the dependence of Lg on because jz is known, E is known, 
and it is not necessary to carry the terms in Lg necesszuy to determine E. 
These simplifications produce 
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The first variation of R, with (2.3.6) as the Lagrangian, with respect to \g 
gives the Euler equation; 


djrAr) _ Q 
dr 


(2.3.7) 


The solution for Ar is simply a constant over r. For Ar to remain bounded 
as r goes to zero, the constant must be zero. The first variation of R with 
respect to Ag and integration by pzirts yields the following Euler equation 
for Ag 



The solution to this equation is A® = Cir + C 2 T~^ where C\ and C2 axe 
constants of integration. Again C2 must be zero for Ag to remain bounded 
as r approaches zero. If Ci is nonzero, a magnetic field will exist in the z 
(longitudinal) direction. For the problems considered here, the only magnetic 
fields that exist are generated by the current I in the wire and Ci is also 
chosen to be zero. 

Because Ar and Ag are identically zero, it is not necessaury to carry 
the terms in (2.3.6) dependent upon Ar and Ag. Consequently, the expres- 
sion for the gauged Lagrangian for the one-dimensional, time-independent 
axisymmetric conductor with a known CTirrent density distribution is 


(2.3.9) 

Notice that for this particular geometry, with time-independent fields, the 
gauge choice for A does not contribute to the Lagrangian and A is completely 
determined by the boxmdary conditions. For this particular case, L is equal 

to R g . 
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The new expression for R is 


R = 




(2.3.10) 


The first ^-axiation of R with respect to Az and integration by parts produces 




dAz 

dr 


(2.3.11) 


where F is the surface of the integration volume considered in the analysis, 
and Tj and rj axe the inner and outer radial limits respectively of the integra- 
tion volume. For this problem, dF is simply dOdz. Substituting the relation 
for B from (2.1.9) {i.e., B = VxA) into the Euler equation in (2.3.11) gives 
the following Maxwell relation and verifys that (2.3.9) is the correct form for 


' 9 - 


I dirBe) 
T dr 


= Mz 


(2.3.12) 


2.4 THE TWO-DIMENSIONAL AXISYMMETRIC CASE. 

The next simplest problem with which to test the four-potential 
method is the two-dimensional zixisymmetric case. As in the one-dimensional 
case, the current is steady (time-independent) and known, p is still zero, and 
cylindrical coordinates are chosen with the rotational axis coinciding with 
the z axis. The four-potential method is now extended to cover this problem 
by eJlowing <)> to vary in the radial and longitudinal directions and Cj but 
not in the circumferential direction e^. Here, and in the sequel, e^, e^, and 
Cz are defined as the unit direction vectors in the r, 6 and z directions respec- 
tively. All partiaJs with respect to 9 now disappear but paxtials with respect 
to z now remain. Since the problem is time-independent, and j is known. 
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partials with respect to t and partials containing $ can be eliminated. The 
gauged Lagrangian is now 


L9 






+ Ag - 


1 d{rAr) dA 


dr 


+ 


dz 


— (jr^r + jeAg + jzAz) 


(2.4.1) 


Note that this Lagrangiain involves all components of A although the inde- 
pendence from d has introduced some simplifications with respect to the full 
three-dimensional case. 

Variation of the above with respect to Ar and integartion by parts 
produces 


<«(«• 


d^J dr 


+ f dT2SAr{r\g^ 

Zi JT2 


(2.4.2) 


where dTi and dT 2 are defined as rdrdO and d9dz in the amd Cr directions 
respectively and Zi and Zj are the lower and upper limits of integration 
respectively, of the integration volume in the Cz direction. To verify that 
the first three terms of the volume integral in (2.4.2) represent a Maxwell 
equation, the expression for B in terms of A$ and Az is needed. The 
correct expression for this problem is 



f > 

Br^r 



B=< 

Bgeg 


(^ - ^) 


Bz^z 

< > 


N 
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The following expression for VxB in terms of Br, B$ and Bz is helpful for 
verification purposes: 


/ 


VxB= { 





(2.4.4) 


Comparison of the 6 component of B with the Euler equation of (2.4.2) 
verifies that it is the Maxwell equation VxB = /xj in the Cr direction. 

Variation of (2.4.1) with respect to Ag and integration by parts pro- 
duces 




df^Ag d 

dz^ "** dr 


dAg \ 
H dz j 


+ 




(2.4.5) 


Comparison of the r and 6 components of B in (2.4.5) verifes that the Euler 
equations match the desired Maxwell equation in the eg direction. 

Finally, variation of (2.4.1) with respect to Az and integration by 


parts produces 


SR{Az)= [ 

Jv 



dV6Az 

dTiSA 



dXg 
■ dz 
dAz ^ 
dr j 


(2.4.6) 


Comparison of the 6 component of B again verifies the derivation of the 
correct Euler equation, this time for the Bz direction. 
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2.5 SUMMARY. 

In this chapter, the terminology and basic backgroimd for dealing 
with EM fields is developed. The four-potentieJ method is also introduced 
and modified for arbitrary materials. The basic Lagrangian for EM field 
problems is presented eind a gauged form of this Lagrangian is also shown. 
To show the broad range of applicability of the fomr-potential method to EM 
field problems, the gauged Langrangian for two simple time-independent 
cases is derived. The first variation of this Lagrcingian, integrated over the 
independent variables, is also taken. This variation is performed to verify 
that the Euler- Lagrange equations for these two particular cases match their 
respective Maxwell equations and zJso to determine the natural boxmdary 
conditions for each case. 

In the next chapter, the four-potential method is extended again to 
obtain the appropriate Lagrangian for two special cases. These cases are a 
conductor with an unknown current density vector j and a conductor in the 
superconducting state. 


CHAPTER III 


CURRENT DENSITY PREDICTING FOUR-POTENTIAL THEORY 

In the previous chapter, the current density distribution j is known. 
Unfortunately, for the general case, neither the path that the current I takes 
through a conductor nor its distribution is known. In this chapter, two 
different cases where I is"known, but j is not, are examined. The first 
case is a normal conductor and the second case is a type I or II (Ginzburg- 
Landau) superconductor. Both cases have an identical geometry, that of a 
one-dimensional infinite wire and both are time-independent with p equal to 
zero. Cylindrical coordinates are used to describe the problem with the z 
cixis coinciding with the rotational axis of the wire. 

The piupose of this chapter is to develop the Lagrangians for each 
of the two problems, and their residuals (Euler equations), so that they may 
be extended to a finite element formulation. Also included in this chapter is 
a brief presentation of the basic theory of superconductivity for types I and 
II superconductors. 

3.1 LINEAR CONDUCTORS. 

The previously derived Lagrangian for the time-independent case in 
three dimensions is 

L, = J dV A)’’ (Vx A) - - f A. + A, (V ■ A)| 

(3.1.1) 



where the superscript T represents the transpose of the matrix or vector. 
The constitutive equation for a linear conducting medium is 
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j = crE = — (tV$ 


(3.1.2) 


As a first guess, (3.1.2) is used to eliminate j from (3.1.1) in terms of the 
variable $. The resulting equation for the Lagrangian is 


Lg = |^(VxA)^(VxA)-ieV$^V$+c7V$^A + Aj (V-A)j 

(3.1.3) 

Integrating Lg over the volume and taking the first variation yields, after 
integration by parts, the following equation 


8R= I dV8k^' 

Jv 

+ J dT^A^ (VxA X n) + (-n Aj)l + J dTd# {eV$ + crA) • n 
" ’■ (3.1.4) 

where n is the unit outward normal to the surface of the volume of integra- 
tion. 


-Vx Vx A + <rV$ — VA, 


■}-/. 


dV8^V • {eV$ - a A} 


The first volume integral is an augmented form of Maxwell’s equa- 
tion VxB = j, whereas the first boundary integral ensures that the B field 
component parallel to the surface is continuous across boimdary surfaces. 
If a is constant across the volume of integration, then the second volume 
integral is a restatement of the Mcixwell equation V • D = 0 because V • (a A) 
= <tV • a = 0. The second boundary integral enforces the condition that 
the normal component of D be continuous across boimdaxies. For the one- 


dimensional problems studied here where the value of e does not change 
across boimdaries, this condition automatically satisfies the homogeneous 


Maxwell equation VxE = 0. 
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If (7 is constant across the whole volume of a conductor, this formula- 
tion presents no difficulties. However, if <r changes as a continuous (smooth) 
function across a conductor, the second Euler equation is incorrect. This 
czin be corrected by augmenting the Lagangian with the constraint V • D = 
0 or by changing the gauge constraint to V • (<7A) = 0. If the conductivity 
changes slowly across the conductor, the conductivity can be approximated 
by a series of step functions. At low temperatures, for the conductors exam- 
ined in this work, the conductivity does change slowly across the conductor 
volume eind the step function approximation is used. This formulation also 
has problems. The second boundary and volume integrals in (3.1.4) combine 
to produce a series of n — 1 equations for n \mknowns where n is the number 
of differing regions that E field passes through. These regions are caused 
by the choice of integration volumes and changing EM material properties. 
Augmenting (3.1.4) by the cmrent conservation constraint, / — /j, dTnc • j, 
where fie is the directed unit normal to the surface that the current flows 
through, solves this problem, fie is aligned in the direction on current flow. 
The new functional Rgcc is 
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Variation with respect to A and $ produces 


8R 


gcc 


= J VxVxA + crV$ — VAjj^ 

- / { V • (eV$ - <tA) + VAe} 

Jv 

+ y dr^A^ |^(VxA X n) + (-nA,)! 

+ J drS^ {(eV$ + <tA) • n + AcnJ 

(^1+ J dr<7nc V$^ 


IT 

+ SXc 


(3.1.6) 


For anticipated extensions to superconductivity, it was originally de- 
sirable to have j as a primary variable whereas the electrical potential was 
of little interest. Rgcc was written in terms of the variables A and j and the 
first variation and integration by parts was performed to give 


SRgcc = J dVdA^ jivxVxA - j - VA,| ~ {ew^j + 

+ J dTdA^ (VxA X n) -h (-nAj)! + Ac ^dTdj • he 


+ 8\c J dThc • 


(3.1.7) 


where u, the resistivity, equeds 1/cr. 

The second Euler equation -f A = 0, which replaces V • D = 
0, is generally incorrect: this is due to the eUmination of V$, which inhibits 
the necessary integration by parts. The lack of this integration also has the 
effect of forfeiting the automatic verification of the homogeneous Maxwell 
equation VxE = 0. These deficiencies can be corrected by augmenting Rgcc 
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with a Lagrangian multiplier field k to produce the new functional, Rp, which 
follows. 


R,=£ dv[-^{VxAf (VxA) - +f^ 


+ Aj (V-A) + k^ (Vxwj 


dThc • j 


(3.1.8) 


At the start of the thesis research, for the finite element work, wj was 
originally substituted for — in (3.1.1). The original equations produced 
a variationcd index of zero for j. This variational index is a constraint that 
was kept as an arbitrary choice to make the research proceed more rapidly 
and results in a formulation that is not the most computationally efficient. 

For the one-dimensional problem, $ as a primary variable, not j, is 
the better choice. A formulation that uses $ is better because it only varies 
in the z direction. This requires only two degrees of freedom over the whole 
domain of the problem to model E and D. With j as the primary variable, 
one degree of freedom per element is needed to evaluate j, zmd a minimum 
of two additional degrees of freedom per element are necessary to evaluate 
K. The j formulation requires three degrees of freedom per dement to model 
the E and D fields. 

Another advantage of the $ based formulation is that with the gauge 
choice V • (cA) = 0, only one constraint has to be augmented to the gauged 
Lagrangian, the current conservation constraint. An additional benefit of 
the $ formulation is that it does not exclude a a that varies smoothly. In 
the thesis formulation, because j is C~^ continuous, a must eJso be C~^ 
continuous to satisfy the homogeneous Maxwell equation VxE = 0. 

However, the formulation that was used to produce numerical results 
here contains j as primary variable and not $ because of time limitations on 
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the thesis research. To reproduce this formulation, it is necessary to integrate 
(3.1*8) by parts to lower the variational index of j to zero. The result is 


By = ^ <iv{ i (V X Af (Vx A) - 
+ A,{V-A)+wj’’(Vx<c)} 

+ j dTujf (« X n) + Ac — y <iTnc • j 


(3,1.9) 


3.1.1 ONE-DIMENfiTONAL LINEAR CONDUCTOR. 

As in section 2.3, the simplest application arises for an infinitely 
long, straight conductor of circtdar cross section. A depiction of the physical 
problem is illustrated in the upper half of Figure 3.1. Again, p equals zero, 
and all partials with respect tp ^ and z vanish. The only nonzero components 
of A and j are Az and jz- By (3.1.2), the only nonzero component of E is 
consequently the only nonvanishing component of k is in the e^ direction. 
The expression for Rp reduces to 




-L 




(3.1.10) 


where dT 2 and dTj are again defined as dOdz and rdrdO respectively. Varia- 
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Figure 3.1: Physical Problem: One-dimensional bulk conductors. 
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tion with respect to Aj, jz Ac and k® and integration by parts produces 



3.2 SUPERCONDUCTIVITY. - ■ - - 

This section presents some of the basic theory of superconductivity 
and the appUcation of the four-potential method to the solution of the time- 
independent superconductor problem. For tHs problem, p is taken as zero, 
and the variable # is no longer required. For cases where E or p are not zero, 
the superconductor behaves as a normal conductor for the E and D fields, 
and these fields can be treated by the methods^ ^scussed in the previous 
chapters. The departure from a normal condjactor is ^chibited in the B 
and H fields and in the resistance of a supercpi^uctor. There is an almost 
complete absence of resistcince zmd the B,^d H fields are non-linear. The 
linear constitutive relation (3.1.2) no longer applies, and j is now a function 
of A and the quantum mechanical quantity, the wave order parameter rp. For 
these reasons, the non-linear fields and non-linear constitutive equations, this 
work deals exclusively with magnetostatic superconductor problems. 

The most widely accepted microscopic theory of low temperature 
superconductivity is due to Bardeen, Cooper and Schreifer [21, pp. 16-71] and 
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is referred to as BCS theory. No attempt is made here to present the BCS 
theory of superconductivity, as the author’s work is based on the Ginzburg- 
Landau equations. The Ginzburg-Landau equations, which describe types I 
and II superconducting phenomena, axe based upon the BCS theory. The 
important result of the BCS theory is that below a certain temperature it 
becomes energetically more favorable for “free” electrons to bind together in 
pairs, called Cooper pairs, and that the density of these pairs in a volume 
can be represented by the quantinn probability density function i/>. Table 
3.1 lists the relevant nomenclature for superconductivity. 

Table 3.1 Superconducting Theory Nomenclature 


Symbol Quantities 

a, /3 Temperature dependent material parameters 

rp Analgous to a wave/position 

function in particle mechanics 
\tp\‘^ Number of superconducting charge carriers 

per unit volume 

Tp* Complex conjugate of tp 

q* Effective charge of charge Ccirriers 

m* Effective mass of charge carriers 

h Planck’s constant divided by 2 tt 

A Magnetic potential vector 

B Total magnetic field 

j Current distribution 

F, Helmholtz free energy of superconducting state 

Fn Helmholtz free energy of normal state 

AF F, - F„ 
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3.2.1 THE HELMHOLTZ FREE ENERGY FOR A SUPERCONDUCTOR. 

The Helmholtz free energy of a system is expressed as 


F = U-TS (3.2.1) 

where F, U, T and S represent the Helmholtz free energy, the potential 
energy, the temperature and the entropy of the system respectively. 

In the general vicinity of the transition or critical temperature for 
a type I or II superconductor, the difference between the Helmholtz free 
energy of the superconducting and normal states of a conductor can be ap- 
proximated as 

AF=F,-F„= f iv{-a|^P + iWI" + 2^lHRV-,-A)VP + iB.H} 

( 3 . 2 . 2 ) 

in S.I. units [20], where the quantities a, ^ and xf) axe defined in Table 1. 
The first two terms represent a typical Landau expansion of the Helmholtz 
free energv for a second order phase transition. The third term represents 
the total momentum of the charge carrier. The —ihV term is analogous to 
the dynamic (kinetic) momentum of a quantum, wave-like particle; the ?*A 
term represents the field momentiim [31, p. 633; 21, pp. 105-108]. 

REMARK 3,2.1 

A good example to illustrate quantum kinetic momentum is provided by a one- 
dimensional particle in an infinitely deep energy well. The — tfiV term in the 
above functional is similar, in quantum theory, to the momentum of the particle 
in the well. 

Using the identities, B = /ioH, and B = VxA, the last term of 
(3.2.2), which represents the field energy, can be replaced by 

(VxA)^ 


1 

2/^0 


(3.2.3) 
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In (3.2.3), the material’s magnetic permeability fi, has been set to 
fio, the value of the permeability of free space. The justification for the 
use of /io is that, in a superconductor, there is an almost total expulsion of 
the magnetic field B from the interior of the superconductor. This effect is 
called the Meissner effect. The B field will only penetrate a small distance 
into the superconductor. This approximate penetration depth is called the 
London penetration depth. For superconducting samples with dimensions 
much larger than the London penetration depth, the contribution to F by 
the difference between and /ioH is small and the substitution of fi for pio 
is justified (Ref. [21], p.89). This type of superconductor is referred to as a 
bulk superconductor. Superconductors with macroscopic dimensions on the 
order of or smaller than the London penetration depth should use fi instead 
of fio‘ Only bulk superconductors are dealt with here. 

Expanding AF in terms of tJ; and ip* gives 

AF=y dv\^-aiPiP* + L^{iPiP*Y + -^{-ihV'^rP-q*A^iP) 

{ihViP* - q*AiP*) + ^(VxA)^(VxA)} (3.2.4) 

Zflo ^ 

The quantities ip and ip* axe both complex quantities and present no 
mathematic difficulties when deriving a variational formulation of supercon- 
ductivity, but they do cause numerical problems. If ip and ip* axe used as 
independent variables in a numerical model, they require twice the amount 
of memory to store because both a real and imaginary number must be 
stored for each variable. A preferred numerical formulation will only contain 
variables that are real. Luckily, the independent variables ip and ip* can 
be expressed in several different manners, all of which are mathematically 
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equivalent. If we let ^ equal tjffi + and xp* equal tpn where tl^ji and 

rf;j represent the magnitudes of the real and imaginary parts respectively of 
the old variables, and i the square root of —1, the memory storage problem 
is solved and the new variables are real. This formulation was used in in 
Ref. [32] for one- dimensional calculations. Although reasonable results for 
most quantities were obtained, others lacked accuracy. Later, it was decided 
to find an improved formvJation. In the modified formulation, ip and ip* 
become l^|e*^ and \ip\e~'^ respectively, where |^| is the modulus and tn is 
the phase angle of ip and ip*. These are the new independent variables used 
in the functional AF. With these substitutions, (3.2.4) becomes 

%/ 

+ — ?* A^) (Wu7 — 9* A)) (3.2.5) 

+ jL(VxAf(VxA)} 

The first variation of AF with respect to \ip\ is 

dy^|V>l{-2a|V'| -1- 2y9|^|* - 

+ — q*A^){JiVzj — q* A))^ (3-2.6) 

+ (<ir%||^nvw| 

The first variation of AF with respect to m is 

« AF(fc) = - dVfo { V . (l^P - 0 a)) } 

+ dTfe |n ■ (l^p - ^a)) I 


(3.2.7) 
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The first vciriation of AF with respect to A is 

SAF(SA) = J dV6A^ 7 (VxVxA) 

— (VxA X n) 
l/^o 

Comparison of the above equations with the Maxwell eqaution V x B = j 
shows that the constitutive relation for a superconductor is 

(3.2.9) 

where j is now a function of A instead of E. Note that j and the constitutive 
relation axe already contained in the Euler equations and that j and Ac are 
not needed as separate variables to mzike the set of equations determinate. 

The set of Euler equations obtained by the variation of AF is collec- 
tively called the Ginzburg -Landau equations. They describe the behavior of 
type I and II superconductors. In the London approximation, V’ is assiimed 
to be constant throughout the conductor volume. For this approximation, 
equations (3.2.6) and (3.2.7) become zero and equation (3.2.8) becomes 

SAF{6A)= f dy6A^(lV»|^^A-b— (VxVxA)| (3.2.10) 
Jv t ^ J 

This type of conductor is known as a London type superconductor. Type I su- 
perconductors are commonly referred to as London superconductors because 
^ is constant over the majority of the conductor voliune and (3.2.10) can be 
used to get a good approximation of the B field inside of the conductor. 

For the Ginzburg-Landau bulk superconductor, rj} becomes a con- 
stant within the superconducting volume at the interior boundary. This 
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means that |V’| is a constant there, and the interior boundary integral of 
(3.2.6) is zero. 

Although the curl of A has been defined, the divergence of A is 
arbitrary. A common choice and the one used here is the London gauge, 
V • A = 0, which is equi'ralent to the time-independent Lorentz gauge. For 
this gauge choice, A must go to zero inside of a bulk superconductor [32, 
p.l2]. tj; must also go to zero at the exterior free space/ conductor boundary. 
This reference shows that, with the London gauge, also be zero at the 
exterior boundary. This condition is equivalent to VjV’l being zero on the 
exterior boundary. With this condition, the outer boimdary integral of(3.2.6) 
is also zero, and the boundary term disappears completely. 

Because of the London gauge choice and the condition that jV’l is 
constant deep in the bulk layer, the Euler equation of (3.2.7) becomes, in the 
bulk region, = 0, requiring that Vm be a constant. The value of the 

constant is determined by energy considerations. The term \xf;\/m*{hVw — 
5 * A) represents the net exchange of field momentum from the magnetic field 
to the kinetic momentum of the charge carriers. Only in the boimdary layer 
is there an exchemge of momentum and in the bulk of a superconductor this 
term must be zero. Because A is zero in the interior of bulk superconductors, 
Vu 7 must also be zero or there will be an exchange of momentum. Therfore, 
for the London gauge choice, w is constant [21, p.l07]. This reduces the 
number of independent variables from three to two. The correct augmented 
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functioned for the generalized three-dimensional case is therefore 



and its first variation is 


6AF,= f <n^S\^\^-2a\^|;\ + 2/3\^f-^V^\^|^\ + \rp\^A' 

Jy m Tn 

+ / -^(VxVxA)-hVA,| 

Jv [ m* /io J 

4- J dT(5A^|i(VxAx n)-l-(-nA5)| 


A^A 


(3.2.12) 

3.2.2 ONE-DIMENSIONAL SUPERCONDUCTORS. 

For the one-dimensional Ginzburg-Landau superconductor that has 
the same geometry as the linear conductor examined earlier in this chapter, 
and no static change density p, (3.2.11) reduces to 


= + 2^ (ft" (^ 


(3.2.13) 





and the first variation is 


6AFg = J dvs\ip\{-2a\tp\ + 2l3\rpf 


An illustration of the physical problem is shown in the lower portion of Figure 


For a London superconductor, ^ is constzint and (3.2.11) and (3.2.12) 


become 


SAF, = -l^SA.{-L(l^(r^-^ 




For both cases, the only nonzero component of j is in the direction 


jz = 

m 


(3.2.16) 
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3.2.3 F.VAUIATTON OF MATERIAL PARAM ETERS a AND B. 

The following is a summary of Tinkham’s derivation of a and /? that 
is presented in [21, pp.105-109]. Appropriate changes have been made to 
convert this derivation from CGS units to SI umts. 

Deep within a superconductor, due to screening. effects (the Meissner 
effect), there are no fields or gradients. The last terms in the functional AF 
drop out and the resulting equation is 

AF = -alV’p + (3-2.17) 

Near the second order phase transition, at the critical temperature 7^, the 
miniTmiTn value for the free energy occurs when 

^ = -2a|^| + = 0 (3.2.18) 

dip 

from which 


\i>? = llfool" = J (3.2.19) 

where \ipoo P is the value for the number density of superconducting charge 
carriers deep within the conductor. Substituting l^ooP back into the preced- 
ing equation for AF, gives 


o?‘ cF’ 


(3.2.20) 


When the critical field Be is applied, AF = — I^Ho. Because of this 
condition, deep within a superconductor, where no gradients are present, 
the following approximation to AF can be made 

2{JLo 2^^ Ho 0 


AF = 


(3.2.21) 
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The work, W, done in setting up a current distribution j [26] is 

W = -\ I f Ady (3.2.22) 

Jv 

From the London theory [21, p. 84], with Ae// equal to the effective London 
penetration depth, the following equation relating j and A can be derived 

1 


J = 




(3.2.23) 


Substitution of this expression for j into the equation for W gives 


W = 


^Iff Iv 


A^AdV 


(3.2.24) 


2/io^e// 

From the Ginzburg-Landau theory [21, p. 107], the expression for the work 
done in setting up a current density j is defined as being 




(3.2.25) 


If gradients of the order parameter are zero 2ind there axe no external fields 
present, the two preceding equations are good approximations to W . Equat- 
ing these two expressions for W gives 

1 2 
2tJLoKff ^ 


(3.2.26) 


Algebraic manipulation produces 


I / |2 _ _ Z 


(3.2.27) 


Solving for ^ ^ves 


^ = OL 


m* 


(3.2.28) 


From before 
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y-o P 

Substitution for P finally yields 


(3.2.30) 

Allowing IV’ooP to equal the number of superconducting electron 
pairs, it is seen that to be consistent with the London theory 

q* = — 2e = twice the electron charge 

m* = 2m = twice the electron mass 


^ ^ r 2 i 2 




(3.2.29) 


3.3 SUMMARY. 

In this chapter, the boundary conditions and the appropriate forms 
of functionals based upon the four-potential method axe determined for two 
conductors with an unknown current density vector. The two types of con- 
ductors considered are a normal linear conductor zmd a superconductor. The 
only approximation made for the hnear conductor is that both u and j be 
step functions. The more general case where they aie both C° continuous is 
also discussed. 

For the superconductor, the Ginzburg- Landau and London type su- 
perconductors are discussed. The London type superconductor is shown to 
be a simplification of the Ginzburg-Landau superconductor based upon the 
assumption that the quantum mechanical variable rp becomes a constant 
throughout the conductor volume. 
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Also determined Eire the boundary conditions for the gauge choice 
V • A = 0. This particular gauge choice reduces the number of independent 
variables for the Ginzburg-Landau superconductor by one. The appropriate 
expressions for the two material parameters for this conductor, a and are 
also determined in this chapter. 

The primary assumption in determining a functioned for the Ginzburg- 
Landau superconductor is that the conductor is near the phase transition 
temperature, %. Fortunately, there is some experimental support that the 
Ginzburg-Landau theory is valid in a much wider range of temperature than 
this narrow range if appropriate values for ^eff -®c used. 

In the next chapter, the thermal dependence of the two fuctionals 
derived here is explored. A functional to predict thermal fields is also pre- 
sented. 


CHAPTER IV 


THERMAL EFFECTS 

In previous chapters thermal effects in conductors have been ignored. 
Thermal effects are quite important in superconductivity because they deter- 
mine whether a conductor remains in the normal or superconducting state. 
Thermal fields also affect the current density distribution in normal linear 
conductors. In order to develop more accurate models of the EM fields, it 
is therefore important that thermal effects be included in numerical models 
of these fields. To accomodate the need to model the thermal fields, this 
chapter presents the functionals for two simple time-independent thermal 
field problems, the heat conduction and heat convection problems. These 
sections summarize material presented in references [24] md [33, pp. 90-92]. 

Typically, the EM materizil properties a, /?, e, fi and u; are temperature- 
dependent. The dependence of e on thermal fields for conductors is mild and 
is not addressed here. The thermal dependence of /i is not discussed ei- 
ther because little experimental data for the test material, extremely pure 
aluminum, could be found. The only available datum found was a room tem- 
perature value[34, p.627], and this value is approximately /i(>. If the value 
of remains within an order of magnitude of /x<, (:.e., ~ 10 fXo), tlie for- 
mulations presented in this work experience no numerical difficulties if the 
correct value for is used. Scaling schemes to improve matrix condition 
numbers and numerical stability axe also presented in the sequel, and they 



52 


call be implemented to cover cases where deviates sigmficantly from fio- 
For the above reasons, the values of fio 3^® substituted for fi and e in 

all numerical formulations presented here. The temperature dependence of 
a, 0 and are discussed later in this chapter. 

4.1 THERMAL FUNCTIONALS. 

For a time-independent three-dimensional system, the functional 

Oj = j dV {^VT'^kVT -qT} - J^dr{Q-nT] (4.1.1) 

can be used to model the heat conduction problem [24, p.2]. Q represents 
the heat flux through the boimdary of the integration volume, T the temper- 
ature, k the thermal conductivity tensor, and q the heat generation rate per 
unit volume. For the time-independent case, all of the above are functions 
only of the spatial coordinates. 

The flrst variation of the above equation with respect to the inde- 
pendent \-ariable T is 

SQi = J dVST {V ■ (kVT) - drSTn • {kVT - Q} (4.1.2) 

For linear conducting media, the heat generation per unit volume 
depends upon the current density j and the resistivity of the material w. 
Both k and u for a material are functions of T, but for the purposes of this 
formulation, they axe treated as functions of the space coordinates. When the 
finite element solution process is discussed, this assumption will be treated in 
a more complete manner. For now, it is assrimed that u; is a function of the 
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space coordinates and the appropriate value of q for the time-independent 
linear conductor is [35, p.ll7] 

g=u;j-j (4.1.3) 

For the case of convection heat transfer, the heat flux across a botmdary may 
be expressed as [24, p.4] 

Q = h(T«,-T) (4.1.4) 

where h is the heat-transfer coefficient tensor, which is only a function of 
the spatial coordinates, and Too is the known free-stream temperature. The 
associated variational functional is 

^v = dTTh • {h (Too - iT) - Q} (4.1.5) 

whose first variation is 

= j dr6Tn-{h{Too-T)-Q} (4.1.6) 



For the one dimensional case, the same geometry 2 is that of previous 
chapters is used. An illustration of the thermal portion of the problem is 
shown in Figure 4.1. Again there is a long cylindrical conductor that extends 
to ±oo in the direction, the longitudinal axis of the conductor coincides 
with that of e^ axis, and the conductor carries a steady current I. Due 
to symmetry, there is no variation of T in the ee and directions. The 
functional 0,^ becomes 



(4.1.7) 
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where k is the thermal conductivity, and like a; is a function of r. Qr is the 
heat flux in the radial direction Cr- The first variation of (4.1.7) is 


= - X (krf) + {.f - q/ 


The one dimensional heat convection functional is 


fiv — f dT2^r |/lconv ^ Qr} 

Jr-t 


(4.1.8) 


(4.1.9) 


ri 


where hconv is the convection coefficient. The first variation is 


<50^, — f dT2^Tr {hconp T) Qr} 

Jt, 


(4.1.10) 


For the Euler equation of (4.1.8), integrating with respect to r once 


gives 


dT / .5 , 

fcr— = - / ujlrdr + C\ 

and integration with respect to r twice gives 


(4.1.11) 



dr +Ci 


/ 




(4.1.12) 


where Cj and C 2 are constants of integration. 

A premise to the analysis of the heat conduction problem presented 
here is that k varies slowly across the domain of integration, i.e., between r, 
and Tj. This premise will be true if the finite elements that eire used in the 
heat conduction analysis can be made small enough to model the temperature 
distribution within the conductor adequately. The only limit on the size of 
the elements is machine accuracy. For the cases where the size of the element 
is smellier than machine accuracy, scaling schemes can be employed to move 
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finite element solutions back into the machine’s range. So, theoretically at 
least, the assumption that k varies slowly across an element’s domain is a 
valid eissumption. This assumption proves to be true for the cases that are 
presented later in this work. If k varies slowly, it can be approximated by a 
linear interpolation across an element. This interpolation is 


k = ki + = ki + (4.1.13) 

where ki and kj are the values of k at the inner and outer boundaries of 
integration respectively. 

In the previous chapter, it is assumed that over elements (between 
boundary limits), that cj and jz are approximated as step functions for a lin- 
ear conductor. That assumption is made again here. With this assumption, 
07 and jz become constants over the range of integration of (4.1.12). Sub- 
stituting (4.1.13) into (4.1.12), and using the above assumption of constant 
current density and resistivity, integration of (4.1.12) provides 


T = 



/ rAr 
\Ak 


kjAr^ 

A&2 


In 




+ C-. 


ki^iki 


rAr 


Ar -f Ak 


+ C2 
(4.1.14) 


where In represents the natural logarithm of the argument. 

If o; is allowed to go to zero, as in a superconductor, (4.1.14) becomes 


For this solution to remain bounded as r,- goes to zero, Ci must be zero. 
Now T is an undetermined constant, C 2 , at rj. This provides an important 
boundeiry condition for any cylindrical heat conduction problem that has 
r,- equal to zero, this boundary condition being that dT/dr\^ equal zero. 
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Referring to (4.1.8), the interior boundary integral, with C 2 substituted for 
T shows that 

Consequently, the heat flux and the temperature gradient at r = 0 vanish if 
the solution of T is to remain bounded. 

Letting rj now equal an arbitrary interior point of the conductor, T 2 , 
the exterior boundary integral of(4.1.8) also disappears for the case of oj equal 
to zero. To find the T distribution, another arbitrary point, r^, between T 2 
and the conductor/free space boundary is chosen. The expressions, (4.1.11) 
and (4.1.15), derived from the heat conduction variational principle, are also 
used again. At T 2 , it is already known that dT /dr is zero. Using (4.1.11), 
it is found that C\ again equals zero and using (4.1.15) determines that T 
again equals an arbitrary constant. This constant must be the same as that 
derived for the case where r varied between zero and in order to satisfy 
the C® continuity of T in the variational functional as well as the boundary 
conditions imposed by the first variation of that functional. Because T 2 and 
T 3 are arbitrary, this requires for the case of zero resistivity that T become 
a single constant over the domain of the conductor. 

This constant is determined by the use of equations (4.1.8) and 
(4.1.10). The former states 

Q, = k^=0 (4.1.17) 

Using this information, (4.1.10) gives 


^ 


(4.1.18) 
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where 7^ and Too represent the surface temperature and the temperature of 
the cooling fluid outside of the convection cooling boimdary layer respec- 
tively. This gives the final result, when w = 0 the temperature distribution 
within a conductor is the constant Too, also that the thermal properties 
k and hconv need not be known. 

For the more general case of a nonzero u;, equations (4.1.8) and 

(4.1.10) are used again to find the temperature distribution T. For this case, 

(4.1.11) must also be used. It is assumed here that is zero and Vj is the 
conductor radius, r^- Equation (4.1.11) then gives 

eYT 

krc-^ = -j ujlrdr (4.1.19) 

where is the conductor’s radius. Combining this result with the boundary 
integral of (4.1.8) gives 

Qr = ~— r ujlrdr (4.1.20) 

re Jo 

Using this result and (4.1.10) gives the following equation 

Ts = — ^ r ujlrdr + Too (4.1.21) 

^conv^c Jo 


At the interior boimdary, is equal to zero and the value of Ci of equation 

(4.1.12) is zero. The value for C 2 can also be determined to be equal to Tg. 
The temperature distribution is now 


T(r) 


-fiu: 


ujlrdr \ dr + Tg 


(4.1.22) 


REMARK 4.1.1 

Strictly speaking, the application of the equation (4.1.15) for T to the supercon- 
ductors presented in this work is only approximate. For these superconductors j 
is a function of A and [V’l, both of which are continuous. This makes j C° 
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continuous, and the integration of the first term in (4.1.11) incorrect because j is 
assumed to be C~^ continuous there. However, because a; is a constant and equal 
to zero for a superconductor, the first term of (4.1.11) disappears, and integration 
of the remaining terms gives the previous result, equation (4.1.15). 


4.2 VALUES FOR. THE THERMAL PARAMETERS k .A.ND K 

The first parameter discussed is the thermzd conductivity A: of a linear 
conductor. The thermal conductivity of a superconductor is not necessary 
for the problems studied in this work and will not be discussed. Reference 
[36] gives a semi-empirical formula for the thermal conductivity of a material. 
This formula is [36, p.6a] 


k = 


1 

a'T ” + 


(4.2.1) 


where 


a 


a 


fm- n) 



(4.2.2) 


For these formulas, k is given in watt cm”^ "I ^ where is in degrees 
Kelvin. The constants m, n, a" and /?' were determined by a curve fit to 
experimental data. The values of these constants for well annealed, 99.9999% 
piire aluminum with a residual resistivity of 0.000593 micro-ohms per cm and 
a critical temperature of 1.196 degrees Kelvin are [36, p.9] 

n = 2.0 a" = 4.8 x 10“® 


m = 2.61 = .0245 


(4.2.3) 


These values are used to determine k for all of the examples presented here. 
The thermal conductivity returned by this formula is accurate to within 3- 
5% of experimental values in the temperature range of zero to fifteen degrees 
Kelvin. 



60 


On a microscopic level, all conductors are composed of a lattice 
structure with tightly bound electrons and protons. Some loosely bound 
electrons also exist and are called “free” electrons. There are two types of 
lattice structure interactions, and they determine, in part, how quickly heat 
may be transported through a conductor. This part is called the lattice 
structure’s contribution to the material’s thermal conductivity or the ther- 
mal conductivity of the lattice. The first of the two lattice interactions is 
that due to quantum lattice vibrations called phonons, that can be treated 
quantum-mechanically as both waves and particles, and collisions between 
these quasi-particles. For this interaction, the thermal conductivity is pro- 
portional to T [31, pp. 115-121], and is represented by the second term in the 
denominator of {4.2.1). The second interaction is due to material imperfec- 
tions, such cis a copper ion in an aluminum lattice structure or imperfections 
in the lattice structure itself, such as dislocations. In this second interaction, 
the transport of both phonons and “free” electrons are being affected by an 
imperfect lattice. The net result is that a particle is being scattered by the 
lattice imperfection. For an essentially pure monocrystaline structure, these 
effects can be neglected. This assumption is made for the above aluminum 
sample for k because of its high pxirity and because it has also been well 
annealed to remove lattice imperfections. 

The “free” electrons provide a third means of energy transport and 
may either trrinsport an electrical current, heat or both heat and a cur- 
rent. The rate of transport is governed predominantly by electron-phonon 
collisions. For a conductor, this is the dominant form of heat transport 
and is called the electronic contribution. This contribution to the thermal 
conductivity is c 2 Llled the electronic thermal conductivity. Electron-electron 
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collisions also occur, but are so infrequent that they may be neglected here. 
For the temperattire range of interest, zero to fifteen degrees Kelvin, theory 
predicts that the electronic thermal conductivity is proportional to 1 ^ [37, 

p.204]. The first term in the denominator of (4.2.1) represents the electronic 
contribution to the thermad conductivity and illustrates the excellent match 
of theory to reality for aluminum. 

The second parameter necessary to the computational analysis of the 
problems posed in this work is the heat convection constant hconv Typically, 
type I and II superconductors axe cooled by liquid helium [38, p.l93]. When 
liquid helium is used, the boundary conditions axe not of simple convection 
coohng, but of combined convection cooling and heat transport by thermal 
conductivity. At the low temperatures necessary to induce superconductivity 
in aluminiim, liquid hehum becomes a two phase fitiid. One part of the fiuid 
behaves normally, and the other part becomes a viscosity free (resistanceless) 
fiuid called a superfluid. Not wishing to model the physics of the superfluid, 
as the focus of the present work is to model thermally coupled superconductor 
behavior, a simple, arbitrary heat convection bovmdaxy term was adopted. 
For this boundary term, it is assumed that the conductor is in a normal 
state, the current density j, the resistivity u), and the thermal conductivity k 
are constants across the whole domain of the conductor, and the difference 
between the surface temperature 7^ and the cooling fluid does not drop below 
one hundreth (.01) of a degree Kelvin. 

The temperature of the cooling fluid Too is known and, together with 
the current /, is one of the two independent loading parameters that are 
varied in the computational analysis of coupled phase-thermeil-EM systems 
presented in later chapters. The maximum values of I and % are used to 
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choose the value of hconv These values represent the state of the system 
where the greatest amount of heat transfer occurs. The above choice of 
state ensures that, by use of the formula for hconv presented below, heat is 
always being removed from the conductor by the cooling fluid and that the 
conductor never cools the flxiid instead. To simplify the determination of 
hconv, aJi overall energy balance approach is used below (e.^., see [33, p.92]. 
In the steady state, or time-independent system, the heat energy produced 
by the conductor must equal the amount of heat energy removed by the 
cooling fluid. For the one- dimensional conductor this is 

2~ f f uj\rdrdz = 2'KTc f hconv {% — ‘Tx>) dz (4.2.4) 

Jz{ Jo Jzi 

For a one-dimensional conductor with constant current density, jz = If Trr c^. 
Substituting this expression for jz into (4.2.4), and using all prior assump- 
tions, produces the following expression for hcOTiV 

hconv = 507r“^rc~^u;J^ (4.2.5) 

For this equation, uj is evaluated at "Ts which equals 7^ -f .01. This choice 
for u> generates the largest possible amount of heat in the conductor for the 
two loading parameters. 

4 3 ' TT^FMAL PROPERTIES OF a;. 

Like the thermal conductivity, two primary mechanisms participate 
to produce a resistance to EM energy transport. For this type of energy 
trEinsport, the electron-phonon interaction predominates again, but only 
dominates at high temperatures (above ~ 20" K). Unlike the thermal prob- 
lem, lattice imperfections can contribute enough to the resistance of EM 
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energy transport that they must be accounted for. The first type of inter- 
action is accounted for by the ideal resistivity of a material. The second 
interaction is accounted for by the residual resistivity of the material. The 
total resistivity is thus the sum of the residual and ideal resistivities 

w = o;,-+a;o (4.3.1) 


where Ui and are the ideal zmd residual resistivities respectively. 

Usually, the residual resistivity is a property of the particular sample 
and is determined by experiment. The value used in the numerical examples 
contained herein is given at the beginning of the previous section. The 
following discussion of ideal resistivity is a summary of material presented 
in Refs. [37] and [39]. 

The ideal resistivity can be expressed as 

ik) ( t ) 


where 77 is a material constant, is the Debye temperature as determined 
by resistance methods, and is 


Js 




dz 


(4.3.3) 


l)(l-e-) 

This is the Bloch-Griineisen formula for the ideal resistivity of a material 
[37, pp. 189-190]. For materials at low temperatures, t.e., 7 r/T ^ 1, the 
upper boimd on J7s can be extended to infinity with little error. Integration 
by parts of (4.3.3) with this new limit produces 


^5 (ar) = 


-t-5 


y.oo ^4 

7o 


dz 


(4.3.4) 
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The first term evaluated at the limits is zero and the second term is 5! x Z (5) 
where Z(n) represents the Riemann zeta function of argument n and 5!Z(5) 
is approximately equal to 124.4. Equation (4.3.2) becomes 

Ui « 124.47^ (4.3.5) 

It bas been observed experimentally, that for low temperatures, u> is propor- 
tional to T® [37, pp. 190-192]. This validates the general behavior of (4.3.5). 
For reasons too lengthy to be discussed here (see Ref. [37], pp. 182-202), 
(4.3.2) and (4.3.5) are only good approximations to the ideal resistivity of 
a material. To bring the formula closer to experimental values, It can be 
replaced by an espression quadratic in T [40, p.470]. Equation (4.3.2) now 
becomes ^ 

Ui = (Co + C,T + C,T^}f^) 

where Co , C\ and C^ are constants determined from experimental data. 

4.3.1 VALUES OF CONSTANTS FOR BLQCH-GRUNEISEN FORMULA. 

The value for Tr is documented as 395° K [39, p.lOO, 37, p.l92]. The 
values for 7^ or Co, Ci, and C2 were not found after an extensive literature 
search. Some constants related to Co, Ci, and C2 were found in Ref. [40]. 
Rather than converting these constants, it was decided to do a curve fit of 
the experimental data in the previously cited reference to determine the de- 
sired constants. The software package Mathematica was implemented using 
the “Fit” option. It was discovered that Co, Ci, eind C2 are not constants 
but pciTcimeters dependent upon the annealing temperature, T4. Curve fits 
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with these parameters assumed to be quadratic functions of 'Ta returned the 
following formulas; 

5Co {7a) = 0.0669523686343453 + 0.00006306135275167563 
- 1.320389349752735 x lO""^ Tj 


5Ci { 7 a ) = -0.001133598163601825 + 0.000006634622902885976X4 
- 7.731210566579611 x 10"^ Tj 

5^2 (Ta) = 0.000003186103199486918 - 1.840858520126625 x 10“® 7a 
+ 2.147449451671961 x 10"^^ Xl 

(4.3.7) 


These empirical formulas agree within 5% when compared with the experi- 
mental data of Ref. [40] over the range of 2.21-273.16® K. For the numerical 
examples presented in later chapters, it is assumed that these values can be 
used down to ~ 0° K with about the same accuracy. This assumption is 
justified because experimental observation shows that in this temperature 
range the residual resistivity is the dominant contribution to the total re- 
sistivity. Fi\^e times the value of each constant is presented in the above 
formulas. This removes a factor of five from the function jTs and simplifies 
of the calculation of The determination of Xs is discussed in the next 
subsection. The value of the annealing temperature used for all numerical 
experiments was 548.16® K. 


4.3.2 NTIMEEICAL APPROXIMATION TO THE INTEGRAL Jk. 

In order to obtain vedid values for Wj, it is aJso necessary to have a 
vahd numerical approximation to J^. Although the numerical results pre- 
sented herein he in the range where 7 ]r/T ^ 1, where the approximation of 
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equation (4.3.2) is valid, in the interest of additional accuracy, v7s is evalu- 
ated numericEilly between its actual limits. The evaluation of jTg between its 
actual limits also ensures valid results for should the solution procedure 
inadvertently step into a range where TrjT is no longer much greater than 
one. 

The range of interest for most of the the author’s applications of the 
Bloch- Griineisen formula lies between absolute zero and about one hundred 
degrees above absolute zero. In this range, numerical approximations to the 
integral in (4.3.3) converges slowly. To improve the rate of convergence, an 
equivalent expression is substituted that is composed of the diiference of two 
integrals. Formally, this expression is 


^5 (x) = r 

Jq 

-i: 


' dz 


(e--l)(l-e--) 
= z^dz 


— -/V 


*dz 


(4.3.8) 


l)(l-e-0 


(e--l)(l 

The first integral, as noted before, is simply 5lZ[5], where Z[n] represents the 
Riemann zeta function of order n. An approximation, good to sixteen deci- 
mal places, as determined by the software package Mathematica for 5!Z[5], 
is 124.4313306172044. This is the value used for the numerical experiments 
contained in this work. Integration of (4.3.8) by parts produces 


5 oo ,oo ^4 

J,(x) = 5!Z[5l + — ^ -5y_ 


(Jr 




OO ^4 


(4.3.9) 


- 1 


-dz 


The integral in the above equation is known as a Debye function. Abramowitz 
and Stegun [41, p.998] give the asymptotic approximation to this integral as: 

= + + + + (4.3.10) 

qz _ I ^ \ n ■n? n* j 
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A thirty term series approximation, was used in the numerical experiments 
presented in this work to evaluate the Debye fimction. This number of terms 
enabled the finite element coding to match the Mathematica software results 
for ui to sixteen decimal places. 

For the author’s numerical experiments, the temperature is evalu- 
ated at the nodal points of each element. This is a consequence of the C° 
continuity of the variational functionals presented earlier in this chapter. The 
problem presented by this formulation is that w for each element is only C ^ 
continuous. To overcome this difficulty, it was decided to evaluate w at each 
node, and calculate the mean of the two returned values. This mean value 
is used eis the resistivity of the element. This ensures that a true mean for u 
over the element is represented. If the mean temperature is used instead, the 
resultant value for u does not represent a true mean because, at low tem- 
peratures, w, is proportional to T®. The mean value that is used assumes 
that u)i varies linearly over an element whereas the second does not. The 
assumption of a linear variation here is consistent with the linear variation 
of all other independent variables of variational functionals presented in this 
work. 

4.4 THERMAL DEPENDENCE OF a. 3. AND 

The thermal behavior of the time-independent superconductor is 
governed by the material parzimeters a and /3. For numerical purposes, it 
was found that it was also necessary to know the thermal behavior of lt/>ooP- 
Equation (3.2.30) of the previous chapter shows that a and yS are both func- 
tions of Xeff and Be- Doss gives the empirical thermal dependence of Be as 



68 


[38, p.65] 

He also gives 


the semi-empirical approximation for Ae// 


A,// (T) = K„ (0) 




(4.4.1) 
as being [38, p.52] 

(4.4.2) 


where A*// (0) and Be (0) are semr-empirical constants that represent the 


effective penetration depth and critical magnetic field when the temperature 
of the system equals zero. For high-purity well-annealed aluminum, Be (0) 
equals 99 gauss [42, p.5] and ^eff (0) Is equal to 500 angstroms [42, p.39]. 
Substitution of (4.4.1) and (4.4.2) into (3.2.30) gives 


a=^Be{0fKff{0f 


m 


13 = 




m’ 


-B,(0)'A,//(0) 


1 - 

1 


n 2 


(4.4.3) 


1 -h {TireY 


Equation (3.2.27) gives the relation that |V><x>|^ equals a/^. Substitution of 
(4.4.3) into this equation gives the thermal dependence of jV’ooP? which is: 


IV’oo (T) 1^ 


tn 


-^e//(0) ^ 




(4.4.4) 


Note that as T approaches 7^, the critical field goes to zero. The 
physical interpretation for this behavior is that any field at Te causes a col- 
lapse of the superconducting phase in a conductor. This corresponds to the 
actual physics of a superconductor. The parameter IV’ooP also goes to zero 
as expected. The parameter Ae// approaches infinity however. The physical 
interpretation of this result is that the penetration of the magnetic field into 
the conductor is complete, again in accordance with physical observation. 
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4.5 SUMMARY. 

In this chapter, the variational functionals that describe two meth- 
ods of heat transfer, conduction and convection, are presented. Important 
assumptions about these functionals are that the thermal conductivity k and 
the heat generation terms q, are functions of the spatial coordinates when 
used with these functionals. The thermal dependence of k and uj are also dis- 
cussed, and appropriate numerical approximations for both parameters are 
also presented. Physical constants necessary to the determination of these 
parameters for the test material used in this work, high punty, well-annealed 
aluminum, are also given. 

An important assumption about the determination of the value of a; 
for finite element analysis is also made. In a previous chapter, it was already 
determined for the specific form of the four-potential formulation chosen for 
mnmerical analysis that w be a step function across an element. In order to 
evaluate w, the temperature T must be known. The thermal functionals of 
this chapter returns two values of T for each element, leaving two choices 
of how to determine an appropriate value of w for each element. The first 
choice is to use the mean value of the two nodcJ temperatures to determine 


the elemental u>: 


(e) . T-W 


= Wo -f- Wi 




The second choice is to evaluate w at both nodal temperatures and use the 
mean of these two w’s for the elemental value of w: 
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The latter approach is used here because it more accurately represents the 
mean value of u for ein element. 

In the discussion of the one-dimensional forms of the heat function- 
als, it is shown that, for a one-dimensional steady state superconductor, the 
temperature T is a constant across the domain of the conductor. This re- 
sult is importeint for two reasons. First, knowledge of the values of k and 
hconv for the superconducting state are not necessary and computational ef- 
fort need not be expended to determine them. Second, since T is constant 
across the domain, no numerical analysis is required to find the temperature 
distribution in the conductor. 

With the completion of this chapter, all the necessary tools have 
been developed for the finite element treatment discussed in Chapters V-XI. 
These chapters show specific^ly how to construct specific elements based on 
the four-potential variational principle and their application to the solution 
of thermal, EM, and quantum phase change problems. 



CHAPTER V 


THE CLEMID FINITE ELEMENT 
The first finite element example of the use of the four-potential 
method for determining EM fields is the simplest. It is the example dis- 
cussed in Section 2.3, an infinitely long, straight conductor of circular cross 
section which carries a known, time-independent, imiform current in the lon- 
gitudinal direction. For comparison purposes, the analytical solutions of 
inside the conductor and in free space are discussed first. 

5.1 ANALYTICAL SOLUTIONS TO THE TEST PROBLEM. 

5.1.1 THE FREE SPACE MAGNETIC FIELD. 

In Cartesian coordinates the radizil component of the magnetic vec- 
tor potential in free space can be calculated from the expression (see, e.^., 
[14,26,27,28,43]) 

= ( 5 - 1 . 1 ) 

47 t jy |r| 

where |r| is the distance between the elemental charge jz dV and the point in 
space at which it is desired to find the field potential. The integral extends 
over the volume containing charges. This expression serves equally well in 
cylindrical coordinates. In fact, the transformation of z components is one 
to one if the center of the coordinate systems coincide. 
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As noted previously, the only non-vanishing component of the cur- 
rent vector is jz dTi where dTi is the elemental cross sectional area of the 
conductor r dr d6 and jz is the current density in the z direction. If d£ repre- 
sents the differential length of the wire, then jzdV = jzdTidi = I<U = Idz 
and |r| = Substitution into (5.1.1) yields 

dz 


( ^ r 


Vr^ + 


( 5 . 1 . 2 ) 


This integral diverges, but this difficulty can be overcome by taking the wire 
to have a ffiaite length 2£, symmetric with respect to the field point, that 
is IzLTge with respect to its diameter. Integrating between —C and +C gives 
the result 


Az{r 


) = ^ r 

J-c 


dz 


■s/r^ + 


47T \ 


+ z 


oi:: 


(5.1.3) 


Expanding this equation in powers of rf C and retaining only first-order terms 
gives 

.4, = -(^)lnr+C, (6.1.4) 

where C\ is an arbitrary constant. For subsequent developments it is con- 
venient to select Cl — hirt, where is the “truncation radius of 

the finite element mesh in the radial direction. Then 

With this normalization, Aj = 0 at r = rj. Taking the curl of A gives the 

B field in cylindrical coordinates: 

( 1 dAg dAff 'v 

r^-Sr 


B = VxA = < 


(Bz] 


< Be 

> = ■< 




oz or 

I \d{rAe) idAr I 

r dr ~r~^ J 



0 ] 

> = < 

dAg 

-df ( 

> 

i 0 J 


(5.1.6) 


It is seen that the only non- vanishing component of the magnetic flux density 


= = g (5.1.7) 

This expression is called the law of Biot-Savart in the EM literature. 


5.1.2 MAGNETIC FIELD WITHIN THE CONDUCTOR. 

Again restricting our consideration to the static case, Maxwell’s 
equations in their integral flux form give 


^ H ■ ds = ^ fjt ^3 ■ ds = j y 


nc<ir 


(5.1.8) 


where C is a contour around the field point traversed counterclockwise with 
an oriented differential arclength ds and flcdT is the oriented surface element 
inside the contour. The term for the electric field disappears in this analysis 
because E = 0. From before, it is known that the right hand side of (5.1.8) 
is equal to the normal component of the current that flows through the 
cross sectional area evaluated by the integral. In the firee space case, this is 
the total current that flows through the conductor. But in the conductor the 
amount of current is a function of the distance r from the center. Again us:.ig 
I to represent the totaJ current carried by the conductor, and Vc the radius 
of the conductor, and assuming an uniform current density jz = I/(7rr^), 
the right hand side of (5.1.8) becomes 


I 


j • ficdT = 




(5.1.9) 


Evaluating the left hand side of the integral and solving for Bg gives: 
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Comparing this equation with (5.1.5), it is seen that if // = /xo then Be is 
continuous at the wire surface r = and hcis the value I{2ttTc)- But if 
fx ^ fio there is a jump (/i — fj.o)I/(2Trrc) in Be. 

The magnetic potential within the conductor is easily computed 
by integrating —Be with respect to r: 


A, 


fil 

47rr^ 


+ Cl 


(5.U1) 


The value of Ci is determined by matching (5.1.5) at r = Tc, since the 
potential must be continuous. The result can be written as 




(5.1.12) 


The preceding expressions (5.1.5), (5. 1.12) for A^ can be verified as being 
correct by substituting them directly into the Euler equation of (2.3.11). 


5.2 FINITE ELEMENT DISCRETIZATION. 


5.2.1 CONSTRUCTING EM FINITE ELEMENTS. 

To deal with this particular axisymmetric problem a two- node “line” 
finite element is sufficient. This provides the C° continuity for Az that the 
variational formulation requires. In the following, individucJ elements and 
element properties are indentified by the superscript (e). The two element 
end nodes are denoted by the subscripts i and j. The magnetic potential Az 
is interpolated over each element as 


Az = NAI') 


(5.2.1) 
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Here the row vector N contains the isoparametric finite element shape func- 
tions for A-. The elements of N are only functions of the isoparametric 
paxcimeter ^ which varies between -1 at node i and +1 at node j. The shape 
functions are 

N = i(l-$ 1 + 0 (5-2.2) 


The shape functions are functions of the spatial variable r and the defining 
relation between r and N is 

r = N{j]|}=Nr<'> (5.2.3) 

contains the nodal values for A. and are only functions of the time t, 

A(') = { } (5.2.4) 

Substitution of these finite element assumptions into the previously derived 
Lagrangian, Equation (2.3.9), and then into Equation (2.2.14), yields the 
variationcd integral as the sum of elemental contributions R = HeR^^\ where 






r _ (.) 

dr dr ‘ 



(5.2.5) 


The nodal values are constant with respect to time because this is a 
steady state problem. Therefore, the integration with respect to t disappears. 
Taking the variation with respect to the element node values of A^*^ gives 




(5.2.6) 


This can be written more simply as 


= p(«) 


(5.2.7) 
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where 



Equation (5.2.7) is purposely written in a notation resembling the stiffness- 
force equations of statics. represents a stiffness matrix derived from 

a potential energy variationed formulation, the nodal displacements and 
represents the external force vector. This form clearly illustrates that so- 
lution and assembly techniques developed for finite element mechanics prob- 
lems can be used to solve four-potential based EM field problems. 


5.2.2 APPLYING BOUNDARY CONDITIONS. 

The finite element mesh is necessarily terminated at a finite size, 
which for this test problem is defined as the truncation radius Tt alluded to 
in Section 5.2. In order to make the boimdary integrals of R vanish, it is 
necessary to look at the boundary integrals of (2.3.11). In the finite element 
formulation, the discretized version of these integrals is 






dr 


(5.2.9) 


where <ir 2 is again dddz, and and ^ represent the nodal values 

for A^ at r = 0 and r = rt respectively and numel is the totzJ number of 
finite elements. Simple observation shows that the first boundary integral 
vanishes at r equal to zero. To make the other term vanish, the nodal value 
for A^ at r equal to rt can be constrained to zero. This is the essential 
boundary condition used for this particular problem. 
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5.3 NUMERICAL EXPERIMENTS. 
5.3.1 TEE FINITE ELEMENT MODEL. 


The test problem consists of a wire conductor of radius Vc trans- 
porting a unit current density. For this problem, the finite element mesh 
is completely defined by specifying the radial node coordinates for each ele- 
ment e as = rn ^ and K the mesh contains Nu,ire elements 

inside the conductor, those elements are numbered e = 1,2, ... N^j^ire 
nodes are numbered n = 1,2, ... iV^.vc + 1 starting from the conductor 
center outwards. The first node (n = 1) is at the conductor center r = 0 
and node n = N^ire + 1 is placed at the conductor boundary r = r^. The 
mesh is then continued with Nfree elements into free space to give a total of 
Nxvire + ^'free + 1 nodes and N^ire + ^ free elements. This type of mesh for 
EM field simulation is unique to four-potential bcised numerical methods. A 
single node is needed at material interfaces to model fields as opposed to the 
double nodes of field based simulations. 

For the calculation of the element stiffness zmd force vectors, the 
material permeability ^ and current density jz are uniform over the element. 
Analytical integration over the element geometry gives 


^ r 1 


/(«) 


-1 


-1 

1 




= I (5.3.1) 


where | -h mean radius of the element and = 

— r-®^ the element radial length. For the test example, is a constant 
inside the conductor whereas outside, is zissumed to be imity. For the 
analyticcil solution of Section 5.1.1, this requires that Ho be replaced by one. 
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The longitudinal curent density is jz = //(Trr^) inside the conductor whereas 
outside jz vanishes. 

The master stiffness matrix and force vector are assembled following 
standard finite element techniques. The only essential boundary condition 
requires setting the nodal potential on the truncation boundary to zero, as 
explained in Section 5.2.2. The modified master equations axe processed by 
a conventional symmetric solver, which provided the value of the magnetic 
potential at the mesh nodes. The magnetic flux density B$, which is constant 
over each element, is recovered in element by element fashion through the 
simple finite element approximation 



dAz 

dr 


£ll A 

dr ^ ~ Z(') 


(5.3.2) 


This value was assigned to the center of each element e for plotting purposes, 
although it is a step fimction due to its C~^ continuity. 


5.3.2 NUMERICAL RESULTS. 

The numerical results shown in Figures 5.1 through 5.6 pertain to a 
unit-radius conductor (tq = 1), with the external mesh truncated at rj = 5. 
The element radial lengths, were kept constant eind equal to .25, which 
corresponds to four internal and sixteen external elements. 

The computed values of the potential Az are compared with the 
anzdytical solutions of Sections 1.1.1 and 1.1.2. As can be seen, the agreement 
between analytical and FE values is excellent. The comparison between 
computed values of the magnetic flux density Bg shows excellent agreement 
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except for the last element near the wire center, at which point the FE 
approximation (5.3.1) loses accuracy. 

Figures 5.1, 5.3, and 5.5 are for the case where fiwire was 10.0, and 
Figures 5.2, 5.4, and 5.6 are for the case in which fiwire was one, that is, the 
same as the space surrounding the wire. Figures 5.1 and 5.2 show computed 
and analytical magnetic potentials. The slope discontinuity at r = 1 in 
Figure 5.1 and the jump in Be in Figtire 5.3 are a consequence of the change 
in permeability fi when crossing the conductor boundary. Figures 5.3 and 5.4 
show the computed and analytical magnetic flux densities. Figures 5.5 and 
5.6 show the computed and analytical magnetic flux densities in free space in 
more detail. Note that Figures 5.5 and 5.6 for r > 1 are identical; this is the 
expected result because as shown in Section 5.1.2, the free space magnetic 
flux fleld depends only upon the current enclosed by a surface integral around 
the wire zmd not on the details of the interior field distribution. 

In summary, this finite element performed very accurately in the 
excimple problem and converged, as expected, to the analytical solution as 
the size of the elements decreased. 
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D Analytical 


Figure 5.3: Be vs. r, fiwire. = lO.O.Values shown on the interface r = 1 with 
dark symbols have been extrapolated from element center values 
to display the jump more accurately; this extrapolation scheme 
has not been used elsewhere. 
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Figure 5.4: Be vs. r, = l-O- 
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Figure 5.5: Figure 5.3 for r > =1, Hwire = 10.0, showing free space Be 

field in more detail. " 
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Figure 5.6: Figure 5.4 for t > =1, = 1-0, showing free space Be field 

in more detail. 
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5.4 SUMMARY. 

In this chapter, the case of a simple one-dimensional infinite wire 
is tested. To perform this test, the linear functional of Section 2.3 is dis- 
cretized using standard FE techniques and. appropriate boimdary conditions 
are determined. When the discretization is complete, it can be seen that the 
governing equations axe of a standard form and present no problems to the 
use of steindard FE solvers for linear systems. 

Analytical solutions for the one-dimensional axisymmetric infinite 
conductor are also derived in this chapter. Presented in this chapter au’e 
graphs that compare results obtained from these analytical solutions and 
from the FE model. The two solutions are in excellent agreement except at 
the center of the conductor thereby validating the use of the four-potential 
method for the determination of EM fields. Most importantly, the four- 
potential method accurately predicts the B field across material interfaces 
without any special boimdary treatment, unlike the conventional field based 
methods. 

In the next chapter, the case of a two-dimensioned problem with 
similar boundary conditions and a known current density is explored. The 
extension of the four-potential method to the two-dimensional case is done 
to offer further proof of the validity of this method for EM field analysis. It 
is also performed in order to show the effect of the Lorentz gauge, as the 
gauge effects disappear in the one-dimensional steady-state example. 




CHAPTER VI 


THE CLEM2D FINITE ELEMENT 
In this chapter, the four-potential FE element for cixisymmetric two- 
dimensional problems is developed. The new elements show the relatively 
eaisy extension of the four-potentizd method through the use of Lagrange 
multiplier adjunction to a broader class of problems. This element was tested 
for two different geometries, a one-dimensional infinite conductor, and a 
cylindrical “can” connected to two infinite feed wires on the top and bottom. 

For both geometries, the current density j is known, and the static 
charge density p is zero. The first geometry is the same as that of Chapter 
IV, and is used to provide a check on the element calculations. The second 
geometry is chosen to allow for a variation of B in more than one direction. 
For this geometry, there is no analytical solution, but the results can be 
examined to determine if they are physically realizable. For this reason, this 
chapter begins with a discussion of the construction of the two-dimensional 
axisymmetric finite element. 

6.1 FINITE ELEMENT DISCRETIZATION. 

In the previous chapter, the ungauged Lagrzingian (2.2.13) is used 
to construct one- dimensional axisymmetric finite elements. In the present 
chapter, the four-potential method is extended to include two-dimensional 
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aodsymmetric problems. In doing so, the basic four-potential does not nec- 
essarily satisfy the gauge condition (2.1.10) o, priori and consequently, the 
gauged form of the four-potential (2.2.15) must be used. 

6.1.1 rONSTRUCTTNC; EM FINITE ELEMENTS. 

For the finite element discretization of the two-dimensional case, 
quadrilateral axisymmetric elements defined by their geometry on the r-z 
plane are constructed. These elements are isoparametric with comer node 
points only. Additional construction details are provided in a later section 
of this chapter. 

In the following, individual elements and element properties are 
again identified by the superscript e in parentheses. The element nodes 

are locally numbered i = 1, n, where n is the number of comer nodes 

(n = 4 for quadrilaterals). The magnetic potential components, Ar, Ae and 
Az are interpolated over each element as 

Ar = Ae = A, = NA^'^ (6.1.1) 

Here the row vector N contains the isoparametric quadrilateral shape func- 
tions, which are only functions of the radial and longitudinal coordinates r 
and z 

N=(iVi(r,z) N 2 {r,z) N 3 {r,z) Ni{r,z)) (6.1.2) 

and column vectors A^*\ A^*\ and A^®^ contain the nodal values of Ar, Ae 

and A I respectively, which are only functions of the independent variable t 
A(‘> = (Ari(t) Ar2(0 Arzit) Ar^(t)) 

A^‘) = (An(0 A, 2 (f) Aezit) Ae,{t)) (6.1.3) 

Al*) = (A,i(t) Az2{t) Azzii) AzM) 
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where, as in the one- dimensional time-independent case, the nodal values 
for and a 1'^ become constants. Substitution of the above as- 

siimptions into the previously derived Lagrangian, Equation (2.4.1), and 
integration over the volume of the element yields the variational integral as 
a sum of element contributions R = where 



and again denotes the volume of the element. Varying the above equa- 
tion with respect to the element node values Aj/'^ produces 



Taking the variation of (6.1.4) with respect to A 


(<) 

e 


gives 


(6.1.5) 




( 6 . 1 . 6 ) 


C'2^ 
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and taking the variation of (6.1.4) with respect to produces 

r. f 1 /sN^aN aN^aN 

«jt(6AW) = ar a7^", 

+ ^9dz 

(6.1.7) 

The variation of (6.1.4) with respect to the last independent variable \g is 

aN . (e)^ 




+ 


dz " 


( 6 . 1 . 8 ) 


To facilitate a more compact formulation, the introduction of the following 
matrix notation is used for the stiffness matrix 


K«<') = 


0 K>‘>a.a. 

0 

K^’\.a. 

symm. 


K(«) 


Ar Aj 

0 


K 


(e) 


Az 

0 


(6.1.9) 


where 


■<“— /.r(?rrs; 

(ei w) (;i <'“')• 

(;s M) 




1 aN^aN 


^(e) ar dr 






(6.1.10) 
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The following vector notation is also introduced to give a more compact 
formulation. 



(6.1.U) 


Using the new notation, it is apparent that the fimte element system can 
be again written in the form of the stiffness-force equations of statics, 
= p(«). Assembling these equations in the usual manner will 
produce the discrete finite element equations of magnetostatics, K^u = p. 


6.1.2 APPLYING BOUNDARY CONDITIONS. 

In Section 5.2.2 it is seen that by constraining A, to be zero at r = 
rt, the boundary integrals for <5Ar vanished. Examination of the boundary 
integrals for SA^ in the two-dimensional case, as shown earlier in Equation 
(2.4.6), show that utilization of the one-dimensional constraint will again allow 
both integrals over <£T 2 to vanish. The physical interpretation of this phe- 
nomena is that at a large enough distance from any axisymmetric conductor, 
the field should always be the same as that of a straight wire independent 
of the conductor geometry. This will occur because at a suflBciently large 
distance, ainy effects, such as the end effects of the “can” of the second test 
example, will decay to zero. 

Symmetry conditions also require that dAz/dr equal zero at r = 0. 
This is most easily achieved by constraining dArjdz to zero at r = 0 because 
dAzfdr = dArfdz there. Constraining Ar at the axis to zero fulfills the 
symmetry requirement. Also constraining Ar to zero at z equal to the upper 
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and lower mesh boundciries will mahe the boundary integrals of Equation 
(2.2.2) disappear. 

The application of these boundary conditions removes the rank de- 
ficiencies of the assembled master stiffness matrix. They are not the only 
bomidary conditions that will work, as examination of Equations (2.4.2), 
(2.4.5) and (2.4.6) show, but they are the easiest to derive, being based upon 
simple physical and mathematical arguments. These are the boundary con- 
ditions that are used for the two-dimensional examples presented herein. 

6.2 NUMERICAL EXPERIMENTS. 

6.2.1 THE FINITE ELEMENT MODEL. 

The finite element formulation described in the previous section has 
been applied to the solution of the two test examples described at the begin- 
ning of this chapter. Both problems are treated with quadrilateral elements. 
Each quadrilateral element has four comer points and one interior node. 
These nodes are defined by their radial and axial positions and At 
each comer i, there are three degrees of freedom, namely Ari, A$i, A^,-. From 
these values, the potential components are interpolated with the standard 
bilinear shape functions, which provide the C° continuity required by the 
variational formulation. The centroidal node carries no physical significance 
and is used solely to provide the extra degree of freedom assigned to the 
Lagrange multiplier Thus each quadrilateral element has 4 x 3 -I- 1 = 

13 degrees of freedom. 

For the calculation of the element stiffness and force vectors, it is 
assumed that the permeability and the ciirrent densities are uniform over 
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each element. The desired stiffness matrix and force vector are calculated 
by numerical quadrature using Gauss formulas. The portion associated with 
the potentials is always evaluated with the 2x2 rule. On the other hand, 
three different schemes were tried on the entries associated with A^: 

Full Integration. The same 2x2 rule as for the potentials is used. 

Selective Integration. A one-point rule is used for and ■ 

Zero Integration. The effect of \g is ignored by omitting the integration of 
the associated terms and placing ones on the diagonal. This numerical device 
effectively forces Aj = 0, and thus “releases” the gauge constraint. 

6.2.2 ASSEMBLY. SOLUTION AND FIELD RECOVERY . 

The master stiffness matrix and force vector are assembled follow- 
ing standard finite element techniques. The boimdary conditions are set as 
explained previously. The modified master equations modified for bound- 
ary condirions are processed by a standard symmetric skyUne solver, which 
provides the value of the potentials at the mesh nodes. 

The physical quantities of interest are not the potentials but the 
magnetic flux density B. This is calculated by discretizing the curl of A. 
Since dA{d9 = 0, the magnetic fields become, after discretization, 


( E ^ 


1 = < 


UJ 

1 5(rN.(e) 


The nodal values for B are obtained by evaluation at the Gauss point followed 
by extrapolation to node locations. The average of these quantities is also 
reported as the centroidal value. As discussed below this value is found to 
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be more accurate than interelement-averaged node values. Consequently the 
centroidal value is used to report results. 

For both test problems, the magnetic permeability = iiwirt. "was 
constant inside the conductor whereas outside it the free-space permeability 
was assumed to be unity. The current densities were assumed 
to be uniformly distributed and consequently were calculated by dividing 
the assumed total current flowing through the conductor by the total cross- 
sectional areas of the conductors. 

6.2.3 PROBLEM 1: A CONDUCTING INFINITE WIRE. 

The first test problem is identical to that reported in the previous 
chapter with a one-dimensional axisymmetric discretization. As shown in 
Figure 2.1, it consists of a wire conductor of radius Tc transporting a total 
current of / = 1 ampere in the z direction. This current was assumed to 
be uniformly distributed over the wire cross section. For this problem one 
layer of quadrilateral elements in the z direction, extending from z = 0 
through z = d, was suflicient; here the distance d was chosen arbitrarily. 
The radial direction is discretized with elements inside the wire and 

Nfree elements outside the wire in free space. The mesh is terminated at a 
“truncation radius” ^ where the potential component Az is arbitrarily 
set to zero. Other boimdary conditions are Ar = 0 on the nodes at r = 0, 
z = 0 and z = d. 

The results obtained with r* = 5rc, N^ire = 4 cind Nfree = 10 for 
the potentials are identical to those reported in the previous chapter, thus 
providing a check on the element calculations. The same results were also ob- 
tained with the three integration schemes noted above for the \g term, which 
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verifies that the Lorentz gauge constraint (2.1.10) is automatically satisfied 
by the finite element shape function for one-dimensional magnetostatic fields. 

The computed magnetic flux density B$ at node points was not as 
accurate as could be expected, generally being too large, especially at r = 0. 
The centroidal values, on the other hand, were considerably more accurate 
as regards matching analytical results. Thus for the second problem field 
values at the element centroids are reported. The extrapolation of B to 
nodal locations is a disadvantage of the four-potential variational approach. 
Field based formulations can compute the value of the B field directly at the 
nodal locations while the four-potential method cannot. 

6.2.4 Pl^OBLEM 2: A CONDUCTING HOLLOW CAN. 

The second test problem, shown in Figures 6.1 and 6.2, brings two- 
dimensional features. It is a hollow conducting cylindrical “can” with infimte 
feed wires connected to the center of its top and bottom faces. These wires 
carry a total current of 7 = 1 ampere in the +z direction; this current 
was assumed to be uniformly distributed over the varying cross sections it 
traversed. For the ends of the can, it was assumed that the current flowed in 
the ±6r direction, forcing js to be zero. For the areas where the feed wires 
join the “can” , and the comers of the “can” , it was assumed that the current 
turned ninety degrees and was uniformly distributed. This assumption is 
unrealistic physically, but warranted for the mesh used in this test problem. 
The mesh choice is discussed below. The wire radius Tc 2 ind the can wall 
thicknesses were assumed to be identical. 

Because of the symmetry of the problem it is sufficient to model only 
the upper half z > 0. The results presented here were obtained by using a 
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25 X 25 element mesh of square elements. Within this mesh, the wire as well 
as the can walls were modeled with only one element across the radius or 
thickness, respectively. 

The regular mesh indeed represents an “overkill” for the free space 
while it is insufficiently refined to capture field distribution details inside and 
near the conducting material. This mesh was actually chosen to conform 
to limitations of the three-dimensional plotting functions of the software 
package Mathematica. 

The problem was irun using full, selective and zero integration 
schemes for the freedoms. The magnetic permeability fifrec i^ free 
space outside the conducting material was chosen as unity. For the conduct- 
ing material two different values for the permeability fi = Hwire were tried; 
1.0 and 10.0; the latter to check whether flux jump conditions were auto- 
matically accommodated by the potential formulation. Selective resffits are 
reported graphically in Figures 6.3 through 6.8. Figures 6.3 and 6.4 show 
the magnitude of Bg for ^tu,ve = = 1 obtained for the full and zero 

order integration schemes, respectively. Figures 6.5 and 6.6 show these re- 
sults in contour plot form. Figures 6.7 and 6.8 correspond to fiteire — 10 and 
show the magnitude of B$ from different viewing points. A discussion of the 
results follows. 

The full integration scheme for A^ performed well outside the con- 
ductor. Results were compared with those of the analytical solution for the 
infinite straight wire (the first test problem) to determine whether they were 
physically reasonable. As r becomes large compared to the can cross dimen- 
sion (towards the outer radial edge of the mesh), the answers agreed. This 
is the expected behavior, because as r goes to oo, the general axisymmetric 
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problem should behave as an infinite straight conductor. As one moved to- 
wards the top of the mesh, the solution again approached that of an infinite 
wire as can be observed in Figures 6.3 through 6.8. This behavior is expected 
because els we move parallel to the wire in the z direction, the effects of the 
current in the can ends should tend to zero and the only fEir-field effects 
shoTild be from the total current. The results for the magnetic field within 
the feed wire were not accurate els it did not vanish for r = 0; this behavior 
was due to the use of only one element across the radius and the fact that 
only centroidal values axe reported as noted above. 

The selective integration scheme gave answers of the same general 
shape ELS the full integration scheme, but they only agreed to one or two 
significant digits; these results are not shown here as they are hard to dis- 
tinguish in plots. The zero integration scheme (which in fact releases the 
Lorentz gauge coupling), gave solutions for the field that were larger than 
expected at the conductor boundary stnd a physically unrealizable field in- 
side of the “can”. This field grows sharply as the can axis is approached, as 
shown in Figures 6.4 Eind 6.6. 
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Figure 6 . 3 : Be vs. r and z for fi^ire = 1- Full integration scheme for Xg. 
Intersections of mesh represent element centroids. 
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Figure 6.4: Bg vs. r a n d z for ^w,ve = 1- Zero integration scheme for A^. 
Intersections of mesh represent element centroids. 
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Figure 6.5: Contour plot of B$ vs. r and z for fiwire — !• Full integration 
scheme for \g. Numbers on axes rep resent the number of elemeiit 
centroids traversed from the center of the “can”. Each element is 
.02 X .02 square. All contours are equally spaced and range from 
minumum to maximum values of the field. 



Figure 6.6: Contour plot of B$ vs. r and z for fi^ire — l-O- Zero integration 
scheme for \g. Numbers on cixes represent the number of element 
centroids traversed from the center of the “can” . Each element is 
.02 X .02 square. All contours are equally spaced and range from 
minumum to maximum values of the field. 
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Figure 6.7: vs. r and z for fixvire ~ 10* Full in.tegra.tion scheme for Xg. 

Intersections of mesh represent element centroids. Note sharp 
field jumps on conductor surfaces. 
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Figure 6.8: The same case as Figure 6.7 shown from a different viewing point 
to emphasize how B$ fails to go to zero as r approaches zero 
because of the coarse conductor discretization. 
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6.3 SUMMARY. 

The results of the CLEM2D finite element show that the fom- 
potential variational principle can be applied to a broeider class of prob- 
lems than the simple one- dimensional axis 3 mimetric conductor. Although no 
analytical solution is available for direct comparison, the physical behavior 
of the numerical results strongly suggest that they are accurate. The only 
point where the results are inaccurate are within the conductor itself. A finer 
mesh grading within the conductor can solve this problem, as results of the 
previous chapter illustrate. 

The only truly unrealistic assumption about the second test prob- 
lem was the assumption of the current density distribution. A physical cur- 
rent will in general not make a ninety degree turn and remain umformly 
distributed. To address this problem, the CUPLE series of elements was 
developed. Given a known current J, these elements can determine the dis- 
tribution of the current density j as well as the B and E fields. These finite 
element models are the subject of the next chapter. 


CHAPTER VII 


THE CUPLEID FINITE ELEMENT 
In this chapter, the four-potential finite element for one-dimensional 
problems with an unknown current density vector is developed and tested 
for two examples. Both examples possess the same circuleir-wire geometry 
shown in Figure 2.1, no static charge density {p = 0), and a known current 
I in the positive z direction. In the first example, all elements have equal 
conductivities. This example gives the same type of fields encountered in 
Chapter V and is used to verify the accuracy of computed solutions. For 
the second example, the element conductivities axe allowed to differ. An 
analytical solution to this problem exists and is compared with the numerical 
solution. 

7.1 ANALYTICAL SOLUTION TO THE TEST PROBLEM. 

7.1.1 MAGNETIC FTFJ.D WITHIN THE CONDUCTOR. 

The Euler equation for 6Xc of Equation (3.1.11) states 

1= [ dTijz (7.1.1) 

JTi 

This is the law of current conservation. For the examples presented here, it 
is cissumed that jx and uj are simple step functions where w is known and jz 
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is unknown. Using these assumptions, (7.1.1) becomes 

numel - 

1= V f (7.1.2) 

where a superscript letter or number in parentheses denotes the element 
number and nunxcl is the total number of elements. The £uler equation for 
Sk 0 of Equation (3.1.11) over each volume disappears because ^ and 
are step functions. If /c is constrained at r = 0 and r = Tc, the set of surface 
integrals for Skq of (3.1.11) do not vanish and produces the following set of 
numel — 1 equations relating the 

= jU+l)J‘+l) (7.1.3) 

Insertion of (7.1.3) into (7.1.2) gives 

numel . . . _1 

Y / ) (7.1.4) 

The above is used to determine this value is then used with (7.1.3) 

to solve for the remainder of the imdetermined 
Equation (5.1.8) states 

jli-^Bds = J i-hcdT (7.1.5) 

For this example, /x is a constant over the volume of the conductor. This 
assumption is discussed at the beginning of Chapter IV . For the one dimen- 
sional case, the contribution for each using (7.1.5) is 


= I j/'V 


(7.1.6) 



where < r < and and represent the inner and outer bound- 
aries of j:^^\ For r = r\^\ is zero. Using the principle of linear 

superposition, the total B field is 

< rf 
(7.1.7) 

The \-alue of is computed by integrating (7.1.6) over r 2 md taking the 
negative of the answer, which is 

= - / B,dr = + Co (T.1.8) 

where Co is again an integration constant. For equal to zero, Co is chosen 
as zero. To ensure the continuity specified for Az by the four-potential 
variational principle, must equal A^^j^ when both are evaluated at 

This requires that for AzP = aJP •, Cq equal — /x 
The value of Co for each region where changes can be evaluated in a 
similar manner. Doing so will give the following expression for A^ 
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7.1.2 THE FREE SPACE MAGNETIC FIELD. 

Using Equations (7.1.5) and (7.1.1), the expression for Be in the free 
space outside of the conductor is 

B, = ^ (7.1.10) 

2wr 

The value for the potentied Az is 

A,=- J Bgdr = -!A]^r + Co (7.1.11) 

Use of Equation (7.1.9) to determine Co gives the following expression for Az 

The above result differs from the previous solution of (5.1.5) by a 
constant. This is not surprising because, in the one-dimensional case, A is 
not vmique and is determined solely by the boimdary conditions. For the 
example of Chapter V, A is constrained to zero at r*. For this example, A 
is constrained to zero at r equal to zero. For a one-dimensional bulk super- 
conductor using the London gauge, A must vanish at r = 0 as discussed in 
Section 3.2.1. Because of this boundary constraint on a superconductor, A 
is also chosen as zero at r = 0 for the one-dimensional current density pre- 
dicting case. This choice is made so that numerical coding that implements 
both elements to model the phase transition of a superconductor will require 
only one set of boundary conditions for A. The consequences of this choice 
axe discussed in the subsection on appljdng boimdary conditions. 


7.2 FINITE ELEMENT DISCRETIZATION. 
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7.2.1 CONSTRUCTING EM FINITE ELEMENTS. 

To deal with this particxilax axisymmetric problem, the two-node 
“line” finite element is again suflBcient. Individual elements and their consti- 
tutive properties are denoted by a superscript (e). The element end nodes 
axe also denoted by the subscripts i and j again. and k$ aire interpolated 
over each element as 

A, = (7.2.1) 


The row vector N contains the isoparametric shape functions for the inter- 
polation of Aj and Kg. The elements of N are only functions of the spatial 
coordinate r. and contain the nodal values for Az and Ke and are 
only functions of the time t, and for the time-independent problem studied 
here, become constants with respect to time. Substitution of these finite 
element assumptions into the previously derived variational functional of 
Equation (3.1.10) gives 








(7.2.2) 


where dT 2 and dTi are again dddz and rdrdO respectively. Variation with 
respect to A^'X \c produces the following expression for the 
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elementzJ stiffness matrix 


K“*'* = 




K 


(«) 

0 

0 


AA 

T 

Aj 


k(‘) 


Aj 


K 


(e) 


K 

K 


JJ 

r 

(e)^ 


(e) 


T 


0 

K(«). 

XV j, 
0 
0 


0 

0 


(7.2.3) 


where 

= / dVf’ 


) 


1 9N^aNl 

i(«) dr dr j 




f=-f d' 

^ JvM 




= - f (7.2.4) 

JvM 

= “/.f 




(7.2.6) 


The above expression for ^ is not complete because it neglects contri- 
butions to from the boundary integral over dT2. The discussion of 

this contribution is deferred to the subsection devoted to the application of 
boundary conditions. Following the notation of previous chapters, u is 
expressed as 


aw 

uW = I 

B 


K 


(7.2.7) 


It is important to note that Ac is a global degree of freedom. This condition 
must be met when the elemented matrices are assembled to form the master 
stiffness matrix. 

Taking the second variation of with respect to the independent 
variables produces the tangent stiffness matrix K, which is identical to \ 
This occurs because is only a function of the radial coordinate r and 
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not the independent variables. This fact is important when non-linear solu- 
tion techniques, such as the corrective Newton-Raphson method, are used for 
the thermally coupled superconducting problem. To use the non-linear solu- 
tion techniques, the tangent stiffness matrix is required, and the equivalency 
of and K means that may be used in the normally conducting 

portion with no modification. 


7.2.2 APPLYING BOUNDARY CONDITIONS. 

As mentioned in the previous subsection, ^ is not complete be- 
cause the boundary integral term over dT 2 is not e\'aluated. The expanded 
form of this term is 




- 27T.ff riV2 

I 


(numef ) 
J 


(7.2.8) 


where H is the height of the element, numel is the total number of elements, 
Ni and N 2 are the shape functions of N for the two-node “line” finite element 
and the superscript terms in parentheses represent m element number again. 
Talcing the variation of (7.2. 8) with respect to the independent variables 
^^(nume/)^ ^( 1 ) j{numei) evaluating at the specified values for r gives 




(7.2.9) 


There remains in a volume term that contains as an independent 
variable. The variation of this term with respect to and produces 


«K<‘> = 0 


(7.2.10) 



Ill 


Evaluation of for a two-node element gives 


) (7.2.11) 


Using this equation, for e = 1 and e = numel is 

{numel (7.2.12) 

_ ^jinumel) j.{numel) ^ J 

Addition of the third and fourth terms of (7.2'.9) to the above gives 
= Sji^^ {2xFo;('> (0 r('> ) 4'^} 

^j(nt,mel)j^(numel) ^(numcl) (7.2.13) 






(nume/) 


o> 


(numc/) 

k\ 


Evaluation of at e = 1 and e = numel and the addition of 

the first two terms of Equation (7.2.9) reproduces the transpose of the above 
results. 

These results have three consequences. The first is that for e = 1 
and e = numel, zeros should be inserted in the appropriate positions of the 
stiffiiess matrix to account for the effects of the boimdary integral over T 2 - 
Performing this operation creates a reuik deficiency in the master stiffness 
matrix. A solution to this problem is to insert a one (1) on the diagonal 
element of K at the appropriate degrees of freedom and to then constrain 
and to zero. This is easily accomplished and causes no ma- 

jor difficulties for finite element analysis. Second, we now have a system 
of numel — 1 equations for tzg. The physical significance of this is that a 
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Lagrangian multiplier is aissigned to each of the boundaries in the conduc- 
tor where w may change, thereby ensuring the verification of the Maxwell 
equation VxE = 0. Again, no diffictdties ensue emd the formulation still 
matches the actual physics of the problem. Third, there are now numel — 1 
equations relating the numel degrees of freedom associated with the 
The latter consequence is the most important because it shows that the con- 
straint Ac is necessary to remove the rank deficiency of the master stiffness 
matrix associated with the 

As mentioned earlier, the boimdary condition on has been 
changed so that the interior node of the conductor is constrained instead 
of the truncation node. The appropriate boundary integral of (3.1.11) is 




(7.2.14) 


As assumed in Section 7.1.1 and at the beginning of Chapter IV, ^ is assumed 
to be constant for the examples of this work. Equation (5.1.8) is then used 
to produce the following result 


H0 = -fx-^ 


dAz 

dr 


I 

27rr 


(7.2.15) 


A minimum amount of algebra and the above relation changes (7.2.14) to 

= -HISAz (7.2.16) 

rt rt 

where the first integral is ageiin allowed to vanish at r = 0. 

In Chapter V, the truncation node at r = is constredned to zero. 
This will produce a reaction force at that node of magnitude H I because 
the imposed constraint is an essential boundary condition. On the other 
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hand, at the degree of freedom associated with there is no boundary 

force because the boundary integral vanishs there. In the case presented 
in Chapter V, the reaction force at the last degree of freedom for was 
not necessary for the analysis of the example problem. This information 
is needed here for the determination of the new boimdary forces. If Az is 
constrained to zero at r equal to zero, a reaction force of -ff I is produced 
at the degree of freedom associated with Az\^^- This situation is analogous 
to changing the end constraint for a one-dimensional bar with a point force 
on the free end from one end to the other. 

For this example, it is necessary to achieve the same loading that 
was exhibited for the example of Chapter V when Az^ ^ was constrained. 
This loading will produce the same B fields but different values for A. Again, 
it is easier to visualize the rational for the above statement by again exam- 
ining the exeimple of a one-dimensional bar again. For a one-dimensional 
bar, this woiold require that the same stresses be produced in the bar for 
the different set of displacements produced by constraining first one end and 
then the other. The validity of this comparison is shown by an exzunination 
of (5.2.1). The expression for is the same as that for a one-dimensional 

FE “bar” element with a linearly varying cross-sectional area. Young’s mod- 
ulus has been replaced by and the cross-sectional area is denoted by 

For the forcing vector we see that jz is a uniformly distributed 
constant loading force. 

For our problem, to maintain the same boundary forces that are 
exhibited when xs constrained, a reaction force of —HI is added 

at this degree of freedom, and another reaction force of 7 is added at the 
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first degree of freedom for Az. The second force is added to cancel the force 
of —HI produced when is constrained. 

The above reaction forces can also be used for a one-dimensional 
superconductor. The forces are the same because variation of the terms in 
AjP associated with Az will produce boimdary integrals that are identical 
to those in Equation (7.2.14). The use of (5.1.8) to determine an analytical 
expression for these integrals will still apply because the integral on the 
right hand side of (5.1.8) requires only the knowledge of the current I within 
the conductor and not its distribution. The only limitation for the correct 
determination of the boundary integral of (7.2.14) is that r > rg. 

Finally, one more reaction force appears firom the variation of 
with respect to Ac- This force has a magnitude of —I and is applied at the 
degree of freedom associated with Ac. Consequently, only three non-zero 
values appear in the global external force vector p. 

7.3 NUMERICAL EXPERIMENTS. 

7.3.1 THE FINITE ELEMENT MODEL. 

The finite element formulation derived in the previous section has 
been applied to two test problems described below. Both problems are 
treated vdth one-dimensional axisymmetric elements. Each of these “hne” 
elements has two end nodes and a common shared glodal node. These nodes 
are defined by their axial positions and Each end node has three 

degrees of freedom. The first degree of freedom corresponds to and 
the third degree of freedom corresponds to From these values, the 

components of the magnetic potential and the Lagrangian multiplier Kg are 
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interpolated with the standard linear shape functions, which provide the C° 
continuity required by the variational formulation. The second degree of 
freedom corresponds to on the interior node while the exterior node has 
no independent variable associated with it on the elemental level and is con- 
sidered “empty” . This second degree of freedom has no physical significance 
and is carried on the interior node so that an extra node per element 
does not have to be injected to account for this independendent variable. 
This scheme is used because it matches the format of the STEP ID finite 
element which carries no injected interior nodes. The use of this scheme 
makes downstream coupling of these elements, when modeling the complete, 
coupled EM-thermal problem, more computationally eflBcient. All entries in 
associated with the “empty” degree of freedom axe assigned the value 
of zero. The common shared global node is injected at the end of the fimte 
element mesh. It carries no physical significance and is used solely to provide 
the extra degree of freedom assigned to A^. Consequently, each element has 
2 X 3 -b 1 = 7 degrees of freedom. 

For the calculation of the element stiffnesses, it is assumed that the 
permeability fi, the resistivity w, the permittivity € and the current density 
jj. are constant over the element. The desired stiffness matrix is calculated 
by numerical quadrature using a two point Gauss rule. As mentioned at the 
end of the preceding section, only three non-zero values appear in p and the 
calculation of is not necessary. 
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7.3.2 APPLYING BOUNDARY CONDITIONS. 

The finite element mesh is necessarily terminated at a finite size. 
For the two test problems, the outer radial end of the mesh is defined eis the 
truncation radius r = r<. The outer radial end of the conductor’s mesh is 
defined as the wire radius Tc. Since current is only carried in the conductor, 
the degrees of freedom for jz between Tc and rt are constrained to zero. 
Similarly, the degrees of freedom for k$ between Tc and r« are also constrained 
to zero. The degrees of freedom corresponding to and are 

also constrained to zero as explained in Section 7.2.2. Az is constrained to 
zero at r equal to zero and H I is injected into p at the degree of freedom 
corresponding to At the degrees of freedom corresponding to 

and Ac, —HI and —I are injected into p. The use of the seven degrees of 
freedom format for each finite element results in a rank deficiency of one 
for the assembled master stiffness equations. This occurs because there axe 
only numel but the elemental degree of freedom format used produces 
numel+1 equations when assembled. The last element only contributes zeros 
to the meister stiffness matrix for the second degree of freedom of the external 
{j) node. To remove the rank deficiency, the second degree of freedom on 
the outer node of the last element is constrained to zero. 

7.3.3 ASSEMBLY. SOLUTION AND FIELD RECOVERY. 

The master stiffness matrix is assembled following standard finite 
element techniques. During the assembly, the elemental entries for and 
^^(nume/) modified as discussed in Section 7.2.2. The extemcil force vec- 
tor is assembled by injecting its three non-zero entries as described in the 
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previous section. The remainder of the boundary conditions are set as ex- 
plained previously. The modified master equations, modified for boimdary 
conditions, are processed by a standard symmetric skyline solver, which pro- 
vides the value of Az and Kg at the mesh nodes, Ac at the injected node, and 
the mean current density over each element. 

As in previous chapters, the physical quantity of interest is the mag- 
netic flux density Bg. The finite element approximation of Equation (5.3.2) 
is used again. However, this time Bg is plotted as a step function to avoid 
the extrapolations necessary to determine the value of Bg at Tc. 

The ability of the potential formvdation to model the discontinutiy in 
the B field at a conductor/free space has already been established in previous 
chapters. For this reason, in both test problems /i and e were set equal to 
one (1) inside the conductor and in the free space siirrounding it. The first 
test problem set all of the u;^®)’s to one, and the second problem set each 
to equal the inverse of the element number (t.e., equaled the element 
number). 

7.3.4 PROBLEM 1: EQUAL CONDUCTIVITIES. 

The first test problem is identical to that reported in Chapter V and 
possesses a one-dimensional axisymmetric geometry. As shown in Figure 2.1, 
it consists of a wire conductor of radius Tc transporting a total current J = 1 
ampere in the positive z direction. The elements were given a unit thickness 
in the z direction. The radial direction is discretized with Nyoire elements 
inside the wire and N free, elements outside the wire in free space. The mesh 
is truncated at a “truncation radius” rt- Boundary conditions were set as 
previously defined. 
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The results obtained with r* = 2rc, Nv,ire = 20, N jrte = 20 for the 
potentials differed from those generated by the previous EM finite elements 
of Chapter V by a constant, as expected. These results eire shown in Figure 
7.1. They illustrate what appears to be an almost exact matching of the 
computed solution to the analytical solution. Analysis of the data values 
shows that at r equal to zero the error is about 33 per cent. The error 
declines rapidly to .2 per cent at and even further to .08 percent at r<. 
This error is attributed to the relatively coarse mesh used for the example 
problem. 

Figure 7.2 shows the results obtained for the computed current den- 
sity. The result obtained is lower than the true value by less than one ten 
thousandth of a percent, thus providing a check on the element calculations. 
Because these results were so close to the exact solution, they were plotted 
as a series of points, rather than a line, so that they could be distinguished 
from the exact solution. 

Figure 7.3 shows the results obtained for the B$ field. To evaluate 
how closely the finite element solution matches the exact solution, it must 
be observed where the analytical solution intersects the tops of the finite 
element “steps” . For an exact matching, the analytical solution will intersect 
the middle of the “step” tops. Although diflficult to see, at r equal to zero, 
to approximately r equal to .1, the exact solution moves right of center on 
the “steps”. This mezms that the computed solution is larger than the exact 
solution. The error in the computed solution ranges from 33 percent at the 
center of the conductor to 4.6 percent at the conductor boundary. Outside 
of the conductor, the error trailed off to .02 per cent. The high error at the 
center of the conductor is due to the relatively coarse mesh discretization used 
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for this problem. Finer meshes were tested cind the finite element solution 
converged to the exact solution as expected. 

7.3.5 PROBLEM 2: DIFFERENT CONDUCTIVITIES. 

This problem’s geometry and the values for r^, r,., Nfree 3^d Ny;ire 
are identical to those used in Problem 1. The only difference is that the 
element conductivity is set to the element number. The values obtained for 
the current densities shown in Figure 7.4 are as accurate eis those obtained 
in the first test problem. 

The computed potential shown in Figure 7.5 displays a behavior that 
is different from that exhibited in the first example. The error ranges from 
a maximum of about 33 per cent at Tc to zero per cent at r equal to zero. 
The error at r* is approximately 13.4 per cent. But the primary quantity of 
interest is Be, not A^. 

The behavior of Be is shown in Figure 7.6 and displays much less 
absolute error than A^. The error at r equal to zero is about 33 per cent, at 
Tc .064 per cent and at r< .016 per cent. The reason for such better results 
for Be is that the rate of change of Az is the quantity of interest, and not its 
magnitude. Referring again to Figure 7.5, it can be seen that the computed 
value for the rate of change of Az appears to be close to the analytical value 
for over half of the range of r. This accoimts for the good values of Be that 
occured for r 5* .2. ^luch of the error that occured in the computation of 
Az can be attributed to the large change in w (= for this example. 

From the first conducting element to the last, there occured a 1900 per cent 
change in the value of u>. Put into this context, the errors that did occur for 
the finite element values of Az su’e reasonable. Several more examples with 
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more slowly varying resistivities were performed to verify that this was the 
source of the error. They axe not presented here because Figure 7.3 illustrates 
what occurs in the limiting case where the conductivity does not vary at all; 
the analytical and finite element solutions converge. Although some error 
remains at the center of the conductor, finer mesh discretizations can be 
used to generate computed solutions that lie within a desired tolerance. 

It is recognized that the error at the center of the conductor for 
both test problems appears large. This is because the error measured is the 
absolute error. Other error estimators are available, but this topic is deferred 
to Section 10.2.3 because of similar errors that occm for both the STEPlD 
and LINTID finite elements. As in the problems studied in this chapter, the 
error measured is the absolute error and appears large. The discussion in 
Section 10.2.3 shows that the error produced by using the STEPlD, LINTID, 
CUPLElD and LET ID finite elements to solve EM and thermal problems is 
within acceptable limits cind that the absolute error alone is not always the 
best measure of a computed solution’s accurzicy. 
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Figure 7.1: vs. r, all elements equal conductivities. 
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Figure 7.2: jz vs. r, all elements equal conductivities. Finite element solution 
plotted as points to highlight accuracy of method in matching 
analytical solution. 
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Figure 7.4; jz vs. r, element conductivity equal to element number. Fi- 
nite element solution plotted as points to highlight accuracy of 
method in matching cuialytical solution. 
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Figure 7.5: Az vs. r, element conductivities equal to element number 
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Figure 7.6: B$ vs. r, element conductivities equal to element number 
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7.4 SUMMARY. 

The results obtained in the previous two problems show that it is 
possible to extend the four-potential formulation to the case where the cur- 
rent density distribution is unknown. This is important since this means 
that it is now possible to solve problems where material and geometric non- 
linearities preclude a linear current distribution. It also means that whereas 
before a knowledge of how the current was distributed within a conductor 
was necessary, with this extension of the four-potential variational principle, 
all that is needed is the total current I through the conductor, its material 
properties fi and u;, and the conductor geometry. 

Having shown the validity of this extension of four-potential theory 
to the prediction of electromagnetic quantities, one is now prepared to con- 
struct a nonlinear conductor, the superconductor. This is the topic of the 
next chapter. 




CHAPTER VIII 


THE SUPERCONDUCTING FINITE ELEMENT 

In this chapter, the four-potential formulations of the Ginzburg- 
Landau and London type superconductors are discussed. Both elements 
use for the example problem the geometry of the one-dimensional infinite 
conductor shown in Figure 2.1. Because the London type superconductor is 
only an approximation to the more exact Ginzburg-Landau equations, only 
the computational results for the Ginzburg-Laundau superconductor are pre- 
sented here. We restrict our consideration to the time-independent (static) 
case. 

For both superconductors, the total current I is known, j is unknown, 
and the static charge density p is zero. Because the four-potential method heis 
shown in the past three chapters that it can easily model the conductor /free 
space boundary discontinuity for Be, only the region within the conductor 
is modeled. The stiffness and tangent stiffness matrices for the Ginzburg- 
Landau superconductor contain the independent variables [V'l and A and 
therefore represent a set. of non-linear equations. A short discussion of non- 
linear solution techniques is included in this chapter as well as a discussion 
on how l^'l and A are scaled to reduce the ill-conditioning of the system of 
nonlinear superconducting finite element equations. 

No analytical solution is available for the chosen problem. However, 
numerical results can be examined to determine if they are physically re- 
alizable. -A.S a second check on the accuracy of the results, the B field as 
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determined by the finite element approximation using A can be compared 
to the B field determined by j of the finite element formulation and Equa- 
tion (7.1.5). Because no analytical solution is available, the first topic to be 
discussed is the construction of the one-dimensional axisymmetric supercon- 
ducting finite element. 

8.1 FINITE ELEMENT DISCRETIZATION. 

As mentioned in the introduction to this chapter, the discretization 
of the Ginzburg-Landau equations results in a set of nonlinear equations. To 
solve these equations, expressions for the residual r, the internal force vector 
f, the external force vector p and the tangent stiffness matrix K are needed. 
For this problem, f and r are determined by taking the first variation of the 
governing functional, p by boundary integrals, and K by taking the second 
variation of the governing functional. The relationship between r, f and p is 

r = f — p (8.1.1) 

In this section, r, f and K are determined eind in the discussion of the 
boundary terms, p is determined. 

8.1.1 CONSTRUCTING EM FINITE ELEMENTS. 

For the finite element discretization of a one-dimensional supercon- 
ductor, the two-node “line” element is again sufficient. Individual elements 
are denoted by the superscript (e). As mentioned in Section 3.2.1, may be 
replaced by /x<, with little loss of accuracy and this substitution is performed 
for the superconducting finite elements derived here. The material parame- 
ters a and /? are dependent upon T. Also mentioned in Section 4.1 is that 
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for the steady state heat conduction problem, T is a constant throughout a 
superconductor. Consequently, the material parameters a, ^ and fip are the 
same for every element and the superscript (e) is omitted for these quanti- 
ties. The element end nodes are denoted by the subscripts i and j. Az and 
1^1 are interpolated over each element as 

Az = ( 8 . 1 . 2 ) 


where the symbols and have been introduced to simplify notation 

later. The row vector N contains the isoparametric shape functions for the 
interpolation of Az and \rl^\. The elements of N are only functions of the 
spatial coordinate r. and contain the nodal values of Az and \ip\ 

respectively, and are time-independent. Substitution of these finite element 
assumptions into the previously derived variational functional of Equation 
(3.2.13) gives 


^ 2m* dr dr ^ 2/z<, * dr dr " 


2m* 


) 


(8.1.3) 

Taking the first variation of AFg^^ gives a set of equations, which collectively 
represent the internal force vector f for each finite element. The portion of 
obtained by varying is 


r(e) 


1^1 




(8.1.4) 
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The portion of the internal force vector associated with is 

I N^NAW 

' JvM \ii^ dr dr • m' 

The total internal force for each element is now 


(8.1.5) 


( 8 . 1 . 6 ) 


This expression applies to a Ginzburg-Landau superconductor. For a London 
type superconductor |V’l^ Is constant and equal to IV’ooP- Consequently, the 
only nonzero portion of is 

Taking the derivative of with respect to the independent vari- 
ables produces the tangent stiffness matrix ^ for each element. The ex- 
pression for a Ginzburg-Landau superconductor is 


= 


where 




rK<*>A,A, 

K<'\.i#r 

(8.1.7) 

f 1 dN'^dN 

|/io dr dr 

m"** J 

(8.1.8) 

yv(«) 1 J 

(8.1.9) 


K(0 


MW ~ 


r f/ /^2 r , A T 

If"’ ^ ^ ar ar 

^ ( 8 . 1 . 10 ) 

For the London superconductor is reduces to K^'^a^a*- 

Examination of and shows that an internal in- 

because both of the independent variables, |^| and 


consistency can appear 
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Az, and their derivatives use the same shape functions. This inconsistency 
can sometimes cause a “locking” problem. For the one-dimensional cases 
considered here, this does not occur and is discussed further in Section 8.2.4. 

For mecheinical elements is the displacement field in the ele- 
ment. In non-linear finite elements, are the visible degrees of freedom. 
The nodal degrees of freedom cannot be solved for directly because the 
internal forces are nonlinear functions of which in turn is a function of 
the “loading parameter” The general technique to handle such nonlineeu’ 
problems is to convert the assembled residual force equations (8.1.0) to in- 
crement ed form by differentiating them with respect to a loading parameter 




df df dv 5p 

IT- = or Kw = q 

ac av ac ac 


( 8 . 1 . 11 ) 


where w is the set of incremental rates and q is the loading vector, w 
cind q represent the rate of change of v and p with respect to a loading 
parameter C- The response v(C) is obtained by numerically integrating the 
above equation in conjunction with Newton-Raphson iteration procedures 
2 is described later in this chapter. The purpose of introducing these new 
quantities here is that they axe necessary for the topic of the next section, 
the application of boundexy conditions. In keeping with the new notation, 
for a Ginzburg-Landau superconductor, and are 


} 




-5C“ 



1 


( 8 . 1 . 12 ) 


For a London superconductor, and are 


v«) = {aW} wM = {^} 


(8.1.13) 
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8.1.2 APPLYING BOUNDARY CONDITIONS. 

The boundary conditions for are addressed first. As discussed 
in the latter half of Section 7.2.2, the discrete boundary terms for the CU- 
PLEID finite element are the same as those of a one-dimensional supercon- 
ductor. The only non-zero values for occur at the degrees of freedom 
corresponding to the first and last nodal displacements of Az- They both 
have a magnitude of H I and differ in the direction of their application. In 
the past, I has been used to represent the total current load. It is now split 
into two distinct parts to give 

I=Io + <:Il (S.l.U) 


where lo represents the initial ciirrent and II the loading current. When 
the loading parameter C is zero, the only load upon the system results from 
the initial current load. When C equals one, by convention the system is 
regarded as being fuUy loaded. Using the new notation, the forcing vector 


Pa, is 


'' 1 ^ 
0 


Pa=(^o + CIl)H{ 


0 


1-1 J 


(8.1.15) 


where 1 and — 1 correspond to the first and last degrees of freedom for Az 
respectively, and the vertical dots represent a continuation of zeros over the 
remaining degrees of freedom for Az- The expression for is 


^Az = 


- 


^ 1 
0 


dC 


= -IlH { 


0 

1-1 J 


(8.1.16) 
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The above expressions for and eo"e valid for any one-dimensional 
conductor where the first degree of freedom of Az is constrained to zero. 

As discussed in Section 3.2.1, VIV*] on the boimdaries is equal to 
zero. Consequently, the boundary terms dependent upon |V’| of (3.2.14) are 
zero and make no contribution to p. The expressions for p|^| and q|^| are 
therefore 

P|^| = 0 q|^| = 0 

2 ind the total external force and loading vectors are 

p = (P^.\ 

I P|*i J I qi*i J 

for a Ginzburg-Landau superconductor. To ensure that there are no su- 
perconducting charge carriers in the free space surroimding the Ginzburg- 
Landau superconductor, \xj}\ is constrained to zero at r equal to Tf.. For a 
London superconductor, p and q reduce to p^^ and q^^^ respectively. 

8.2 NUMERICAL EXPERIMENTS. 

The finite element formulation for a Ginzburg-Landau superconduc- 
tor htis been applied to the solution of a one-dimensional axisymmetric infi- 
nite wire. Each element contains two end nodes and a common global node 
thai; is located at the truncation radius r^. These nodes are defined by their 
radial positions and The glodal node carries an “empty” degree of 
freedom and is used only to provide the same number of degrees of freedom as 
contained in the CUPLEID finite element. Similarly, each end node contains 
three degrees of freedom. The first and third degree of freedoms carry the 
nodal values for |^| and Az respectively. The second degree of freedom is also 


(8.1.17) 


(8.1.18) 
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"empty”. This choice provides for easier downstream couphng of the super 
and normal conducting finite elements by reducing computational effort. 

The computational effort is reduced because nodal connectivity and 
freedom tables used to generate the diagonal location pointer array for the 
skyhne symmetric stored system of equations are only generated once. The 
“empty” degree of freedom on each end node also allows for the easy addition 
of the variable zo if the gauge choice used is not the London gauge. With the 
“empty” degrees of freedom, each element carries seven degrees of freedom 
like the CUPLElD finite element. The actual ntimber of degrees of freedom 
used per element is 2 x 2 = 4 degrees of freedom. 

For the calculation of K, p, q and f the permeability n is set to fio, 
8is discussed in the introduction to Chapter IV. The values for a, P and jt/^ooP 
for each element are determined by using the formulas presented in Section 
4.4. The tangent stiffness matrix and internal force vector are calculated by 
numerical quadrature using a two point Gauss formula. 

8.2.1 APPLYING BOUNDARY CONDITIONS. 

The finite element mesh is terminated at rj, as in the linear conduc- 
tor. To ensure that no superconducting flux Ccin cross the conductor’s outer 
edge into free space, the degrees of freedom corresponding to j^| between Tc 
and Vi are set to zero. At r = r^, |^| is also set to zero. By doing this, the 
botmdary terms of (3.2.14) vanish. At r = 0, A^ is set to zero as required 
by the London gauge choice. Any “empty” degrees of freedom are also con- 
strained to zero to prevent rank deficiencies of the assembled master stiffness 
matrix. 
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8.2.2 ASSEMBLY AND SOLUTION. 

The tangent stiffness and internal force vector are assembled follow- 
ing standaurd finite element techniques. The tangent stiffness K is stored 
using a symmetric skyline storage scheme, and then modified for boundary 
conditions. The external force and loading vectors are inititalized to zero 
and the two non-zero values for each are injected at the appropriate degrees 
of freedom. 

Solution Technique. 

For linear finite elements, the displacements v can be solved for di- 
rectly by inverting the stiffness matrix and multipl 3 dng it by the external 
force vector p, i.e., 

v=(k")~'p (8.2.1) 

The standard technique shown above to solve for v does not work for the 
Ginzburg-Landau superconducting finite element because is a fimction 
of IV’I and Az- To begin our discussion of nonlinear solution techniques, the 
residual equations are rewritten as 

f - p = r = 0 (8.2.2) 

where f and r represent the internal force 2 ind residual vectors respectively. 
It can be seen that when the residual vector is zero, the solution vector v lies 
upon an equilibrium curve or path called the response. The central idea of 
non-linear solution techniques is to find a solution that lies upon a physically 
correct equilibrium path and to then advance the solution zilong it. For 
the cases examined in this work, the position along the equilibrium path is 
determined by the loading parameter, also known as the control parameter, 
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The displacements v are also called the state variables because they represent 
the state of the system along an equilibrium path. For the cases presented 
here, an initial solution that lies upon the physically correct equilibrium path 
is not always known in advance. For these cases, we are required to guess 
a “neighboring” state from which to start an iterative process that takes us 
to the path. The initial solution, or “guess”, is named the Teference state of 
the system. 

For the Ginzburg-Landau superconducting finite element (STEPlD), 
it was found that the best choice for the reference state is that where Ai is set 
to 0 and all unconstrained values of |V’| are set to |^oo|- The value for |V>ool 
is determined from the formula presented in Section 4.4. This state closely 
approximates a Ginzbmrg-Landau superconductor with the total current I 
and external B fields equal to zero, the difference for our choice occuring 
primarily in the boimdary layer. The true state can be closely approximated 
by a step fimction for |^| over the interior of the conductor with a magnitude 
of 1^00 1- However, the chosen reference state is close enough to the true state 
that the same techniques used to advance the solution can also be used to 
bring this reference state onto the desired equilibrium path. 

To find how the solution vector v changes as the solution advances 
eilong an eqxiilibrium path, the partial of r with respect to the loading pa- 
rameter C is taken to give 

dr _ df dp 

_ di dv dp 

dv d( d( 



= Kw — q = 0 


(8.2.3) 
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The last equation shown above is known as the incTcmenial Tate equation. 
K defines the tangent to the equilibrium path zind is known eis the tangent 
stiffness matrix; w represents the rate of change of the solution vector along 
the equilibrium path and is named the incTemental velocity vector; and q is 
the loading vector that represents the rate of change of the system’s external 
forces as the external loads of the system are varied. 

To advance the solution along an equilibrium path, the values of the 
solution vector at the current known state axe used to determine the tangent 
stiffness matrix. The loading vector is also determined and the following 
system of equations is solved to determine the incremental rates. 

w = K‘^q (8.2.4) 

Numerical problems arise, however, if the current position of the solution 
on the equilibrimn path is a stationary (critical) point. At these points, K 
is singular. For the STEP ID finite element, this occurs when AFg is zero. 
There is no difference between the Helmholtz free energies for the supercon- 
ducting and normal states of a conductor at this point. This point represents 
a crossing of the equilibrium paths for the Helmholtz free energy functionals 
of the normal and superconducting states and is called a bifurcation point. 
If this point is reached or exceeded, because the free energies are equal, the 
LETID finite element is used. The LETID formulation is not singular at 
this point because it does not contain the qucmtmn parameter |^|; conse- 
quently it does not model the true state of the system at this point. The 
tn?e state is a mixture of normal and superconducting phases that lies be- 
yond the modehng abihty of one-dimensional finite elements as it is in fact 
a multidimensional problem. 
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Assuming that a current value of v on the equilibrium path is known, 
as a first step to obtaining another solution, an increment along the tangent 
to the equilibrium path is taken. That is, the solution is moved Av in the 
V direction and in the ^ direction of the hyperspace defined by v and 
The step is named the predictor step. New values for v and ^ are computed 
at the point that lies at the end of the predictor step. A corrector procedure, 
called the corrector step, is then invoked to iterate the solution back onto 
the equilibrium path. The distances traversed in the v and C directions for 
each iteration are designated as d* and 77 * respectively, where superscript 
k designates the iteration number. The equilibrium path is reached when 
r = 0. To ensure at each step n that the solution does not travel too far 
from the equilibrium path, a distance is specified. The distance is 
also used to ensure that the distance along the equilibrixim path traversed 
is not too large. This distance is limited so that the solution procedure 
does not accidentally step over a stationary point or move too far from the 
equilibrium path. Detecting stationary points becomes important when they 
are branching or bifurcation points because it is desired that the solution 
procedure follow the equilibrium path that matches the true physics of the 
system. If the solution procedure steps over one of these points, it may follow 
a non-physical equilibrium path. 

The addition of the length constraint adds an extra equation to the 
original system of equations: 

|As„|-/n = c = 0 (8.2.5) 

where !As„j is an approximation to the distance s travelled on the equilib- 
rium path. For the finite elements of this work, the initial values for v and 
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^ are computed using a Forward Etiler scheme. The corrector step uses the 
Conventional Newton-Raphson (CNR) method to iterate to a solution under 
the arclength constraint (8.2.5). The formulas used for the forward Euler 
integration and the arclength constraint are reproduced below. Subscripts n 
represent the step number and superscripts k represent the iteration number. 

Arclength Constraint 

lASnl - = -^ |w^Av„ + ACn| - /n = 0 

Jn 

fn = + wT"w„ (8.2.6) 

dc j, w„ dc _ 1 

dv ^ fn dC ^ fn 

Foward Euler Method with Arclength Constraint 
Cn+l “ Cn ”1“ ^Cn ^Cn ^ ^n/ fn 

Vn+I = v„ + Av„ Av„ = K“^q„ACn (8.2.7) 

= w„AC„ 

To implement the Conventional Newton-Raphson technique, the 
original equations for r must be augmented by the constraint equation c 
and solved. This gives the linear system 

[•■5 71W-H 

Because this augmented system is not symmetric, the two lineeir symmetric 
systems below are solved for instead 


Kdr = -r Kd, = q 


(8.2.9) 
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to get 

d = d. + ,d, (8.2.10) 

which finally gives 

v;+‘=vs+d ci« = c„‘+’) (8.2.11) 

Iterations are performed until the 2-norm of r is less than a specified toler- 
ance T. For cases where the 2-norm will not go below t, a limit is set on 
the maximum number of iterations by another input parameter. Because the 
Newton-Raphson technique can also diverge instead of converge upon a solu- 
tion, limits on the maximum value for the 2-norm of v are also specified. To 
stop the solution process, another input parameter limits the maximum value 
of C- ^Vhen this value for C is surpassed, the solution process is terminated. 
The solution procedure may be summarized as follows: 

(1) Initialize v and C to the reference state. 

(2) Solve for w. 

(3) Update v and C using the Forward Euler integration scheme 

(4) E'caluate r and c at step n -t- 1 by using v„+i and Cn-n- 

(5) Solve (8.2.9) for dr and dg. 

(6) Using (8.2.10) with the values for dr and dg, solve for rj eind d. 

(7) Update v and C using (8.2.11). 

(8) Find the 2-norm of v 

(9) If the maximum value for the 2-norm of v or ^ is exceeded, terminate 
the solution procedure. 

(10) Find the 2-norm of r and c. 
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(11) If the 2-norm of (10) is less than or equal to r, restart the solution 
procedure over at step (2) until the desired value for ^ is reached or 
exceeded. If the 2-norm of (10) is greater than r , go back to step (5) 
and repeat the solution procedure until the 2-norm is less than or equal 
to r or the maximum number of iterations is surpassed. 

8.2.3 SCALING TECHNIQUES. 

The solution procedure of the previous section is particularly sensi- 
tive to heterogeneous physical dimensions in the solution vector v. In ad- 
dition, off-diagonal terms of K may be either considerably larger or smaller 
than the diagonal terms, giving K a high condition number. This means 
that a small change in one degree of freedom may produce a large, non- 
physical displacement at another degree of freedom. The STEPlD finite 
element solves the above problems by implementing several different scaling 
techniques. The first technique gives the elements of v the same physical 
dimensions. The second seeding makes off-diagonal terms of the same order 
of magnitude as the diagonal terms. The third scaling is used to further 
improve the condition number of K. Finedly, the fourth scaling adjusts the 
dimensions of v and C in the solution hyperspace to improve convergence 
rates and accuracy. 

To perform the first, third and fourth scalings, the solution vector v 
is scaled by a diagonal matrix S„ 

v„ = S„v (8.2.12) 

where subscript n indicates either the first, third or fourth scaling and the 
superposed tilde denotes a scaled quantity. If the stiffness-force equation is 
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premultiplied by S“^ and the scaled form of v is substituted, the results are 

s;'k"s;'v = s;'p 


K^v = p 


(8.2.13) 


where eind p aire the seeded versions of and p. The scaled versions of 
other relevant quantities can be derived in a similar manner and are presented 
below. 

f = s;'r f = s;'f q = s;'q 
w = s;V Av = s;'Av K = s;'Ks;^ 

(8.2.14) 

|As| = jw^S^Av + AC| = |w^Av + AC| // 


/ = yjl + w^S^w = Vx±wJw 

First Scaling. 

The first scaling is performed element by element at an element level, 
and is used to scale \ii?\ and to have the same dimensions. Let L, M, 
T and Q represent umts of distance, mass, time and charge respectively. 
It can be seen that has units of MLfTQ and 1V»1 has umts of 
Numerical experiments showed that letting the units of v be L ^ improved 
the stability of the solution process. For this scaling. Si for each element is 
expressed as 


o(e) 

*^11 


0 

1 



0 0 0 0 O' 

0 0 0 0 0 

5^3^ 0 0 0 0 

5^4^ 0 0 0 

10 0 

*^66 

1. 


(8.2.15) 


symm. 
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c(®) c(*) 

*^11 ” *^44 


= V^ 


c(*) _ c(®) ^ 

■^33 ~ *^66 ^ 


where ones have been placed at the degrees of freedom corresponding to the 
“empty” degrees of freedom. This allows for the inversion of S[^\ Both 
scaling factors are constants over the domain of the superconductor and do 
not affect the assembled scaled master stiffness equations if this scaling is 


performed at the elemental level. 

Second Scaling, 

The second scaling is performed to mahe off diagonal elements of K 
of the same order of magnitude as the diagonal elements. It also serves the 
dual purpose of bringing steps in the v-^ hyperspace into a more reasonable 
range. The second scaling is essentially a conversion of umts from one system 
to another. After performing the first scaling, the umts of v are where 
L is meELSured in meters, the appropriate umt of length for the rationalized 
MKS system of measurement. However, most of the material parameters 
for a superconductor are of the order of about 10 ® meters. To m a k e the 
order of magnitude of the off diagonal terms approximately the same as the 
diagonal terms, it was observed that on an element level this could be done 
by changing the units of length to micrometers (10“^ meters). To perform 
this conversion, etU of the nodal positions, and are multiplied by 10®. 
The permeability of free space has units of M hj and is also multiplied by 
10®, whereas H has units of M jT and is multiplied by 10^^. The effective 
penetration depth Ae// has units of L sind is also multiplied by 10®. These 
are the only quantities that are changed to perform the units conversion. The 
remaining material parameters a, ^ and |t/»ooP are calctdated using the new 
values for /Xq and Ag// while the scaled value of Ti is used for calciilating Si 
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and B has units oi M/TQ and is not affected by this unit conversion, 

so scaled quantities can be used “as is” for field recovery. 

Third Scaling. 

The third scaling is performed on the assembled tangent stiffness 
matrix K. It is a simple diagonal scaling where nonzero elements of the 
diagonal matrix S3 are equal to the square root of the absolute value of the 
corresponding diagonal element of K, i.e., Sa = y/\Kii\. For the “empty” 
degrees of freedom, the diagonal elements of S3 are set to one to give full rank 
to the matrix so that it can be inverted. This is a common scaling technique 
that will reduce the condition number of a symmetric positive definite matrix. 
Although the constrained stiffness matrix for the superconductivity problem 
is negative definite, this technique works well here. 

Fourth Scaling. 

The fourth scaling is also performed on a global level. Again a diag- 
oneil matrix S4 is used, but this time all of the elements are the same. The 
purpose of this scaling is to make AC approximately equal to I„. To meet this 
requirement, / must be approximately one. It is ensured with this require- 
ment that no matter what the value of the product w^w may be, the scaled 
distance traversed along the equilibritim path will be approximately equal 
to the desired input distance The value for each element of S4 that gave 
the best resiolts for the STEP ID finite element numerical examples presented 
here was 10^. 
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8.2.4 MESH GENERATION. 

The superconductivity problem exhibits boundary layer characteris- 
tics because most of the physics occur in a relatively narrow region close to 
the conductor/free space boundary. The finite element mesh must have a 
fine grading in this region to model superconductivity accurately. Eighteen 
months of research and numerous niimerical experiments have shown that 
if the mesh grading there is inappropriate, the solution generated will suf- 
fer accordingly. The problem most often encoimtered by poor mesh choices 
was that of a high condition number for K. This generally causes the so- 
lution method to fail because the 2-norm of the residual vector r and the 
arclength constraint c cannot be brought below a reasonable value for the 
input tolerance r. The solution “dances” aroimd the v - C hyperspace until 
the maximum number of corrector iterations is reached or the 2-norm of v is 
exceeded. In a few rare c^es, with a poor mesh choice, the solution actually 
did converge. These solutions were rarely of any value because the condition 
number was estimated to be in the to range of 10® to 10^®! Another diffi- 
culty encountered with a poor mesh choice is that |^| and Az will oscillate 
around their equilibrium values. To solve these problems, it is necessary to 
reexamine the theory of superconductivity. 

In the previous discussion of superconductivity, the effective London 
penetration depth Ae// was introduced. The London penetration depth pro- 
vides a measure of how fax the B field penetrates into a superconductor from 
the conductor/free space interface. This is significant because Ag// provides 
a minimum depth for the boundary layer that is being modeled. This is the 


144 


range over which Az should decay to approximately zero. It is also neces- 
sary to know the rzmge in which the other variable, |V’l, decays from its bulk 
layer value of iV’ool to zero. To accomplish this goal, the Ginzburg- Landau 
equations must be examined once again. The following is an abstraction of 
material presented in Reference[ 21], pp. 111-114, and is used to determine 
the range of decay of |i/’|. 

The variation of AFg for a one-dimensional superconductor in Carte- 
sian coordinate space will produce the Euler equation 




(8.2.16) 


where x is the one-dimensionzil spatial coordinate and A has been set equal 
to zero because we axe primarily interested in the behavior of ]^|. If the 
normalized wave function which equals |V’l/lV’oo| is introduced, and 

some algebra is performed, Equation (8.2.16) becomes 


,1/1 I /|3 n 

2m-a +Wn-W/v = 0 


(8.2.17) 


Linearizing this equation by setting equal to 1 -h 6(x), where b(x) C 1, 
gives the first order expansion of this equation as being 

ft" a" 4 (i) ^ _ (1 + (1 + 

(8.2.18) 


2m* a dx^ 
dx^ 


The first term of the equation shows that the decay of is determined 
by I2m* a. This length is referred to as the Ginzburg-Landau coherence 

length i{F). Appropriate substitutions from Chapter IV will give 




,21 


1 + (T/Te) 

1-(T/T,)2J 


(8.2.19) 
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To avoid confusion with the isoparametric coordinate this length shall 
edways be referred to as The dimensionless Ginzburg-Landau param- 

eter k(T) is also introduced, which is the ratio of the two lengths ^(T) and 


k{T) = 


((T) h 





( 8 . 2 . 20 ) 


where «(T) shall be used for this ratio to distinguish it from the Lagrainge 
multiplier vector k. A superconductor with k{T) < l/y/2 is called a type I 
superconductor, while a superconductor with k(T’) >1/ is called a type 
II superconductor. Figure 8.1 shows the difference between and ^eff 
for type I and II superconductors. 

For the particular case of high purity aluminum, k( 0) « .1. This 
maJces it a type I superconductor and shows that the decay depth for |^| is 
approximately ten times the decay depth of Az, where Ae// is the approx- 
imate decay depth for Az- Consequently, the boimdary layer region to be 
modeled must have a depth of at least 10 x Ae// to capture |V’|, furthermore 
^(T) determines the size of the boimdary layer mesh. For a type II super- 
conductor Ae//(T) > ^(T) and the size of the boimdary mesh is determined 
by Ae//(T). Numerical experiments confirm that for aluminum the mesh 
choice of 10 x Ae// reduces the condition number of the system. Numerical 
experiments also show that the mesh generated must be a function of T 
because Ae// and ^(T) are both functions of T. The results obtained with 
the above mesh show realistic values for \rj)\, but both |t/>] and A. exhibit 
oscillatory behavior. Expanding the boimdary layer depth to 200 x Ae// 
caused the oscillations in \ip\ to disappear. All elements generated in the 



Type n Superconductor 

Figure 8.1: Differences between B, ^(T) and A^// for type I and II super 
conductors. 
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boundary layer were equal length elements, where the element length was 
equal to the depth of the boundary layer divided by the number of boundary 
layer elements Nbound- The elements used to model the bulk layer were also 
“regulzir”, their length being equal to (rc — 200 x Xeff{'^))/^bulk, with Nt,uik 
representing the number of bulk layer elements. 

The oscillations are triggered by three different error sources. The 
first one comes from the approximation that is made for (8.2.18). The 
Ginzburg-Landau equations are linearized there to get an idea of the pen- 
etration depth of the magnetic field. The coherence length is only a rough 
approximation to the true penetration depth and not an exact one because 
only the linearized system of equations has been solved and not the exact 
system. The second source of error arises because the finite element model 
is not exact. It merely tries to approximate the continuous case by discretiz- 
ing the region of interest. The third source of error is that fimte precision 
mathematics are used when a solution to the discretized superconductor is 
attempted. The solution procedure and the scaling procedures used all in- 
troduce ntimerical error into the computed solution because of the machine’s 
inability to resolve numbers beyond 16 significant numbers. The expansion 
of the boundary layer helps to push the oscillations induced by the numerical 
error of the solution and scaling techniques below machine limits and more 
importantly, accurately captures the physics of the problem. 

After the oscillations in |^| axe removed, Az may still exhibit oscilla- 
tions close to the conductor /free space botmdary. It was thought that they 
were induced by the mesh being too coarse for that region. Numerical exper- 
iments showed that this was indeed the case. There are three methods that 
can be used to resolve this problem. The first method consists of generating 
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another boundary layer of regular elements neax the conductor /free space 
boundary. The second method involves changing the length of each element 
so that the mesh is more finely graded at the conductor/free space boundary 
than at the interior edge of the bormdary layer. The third method is to 
simply insert more elements into the boundary layer. All three methods add 
extra node points at the conductor/free space boimdzury emd serve to make 
the finite element approximation more accurate. The third method was used 
for the examples here to expedite the research. This is the least compu- 
tationally efficient of the three, but time constraints limited the author to 
using this choice. 

It is mentioned earlier in Section 8.1.1 that the use of the same shape 
fimetions for the calculation of |V’| and A^, and their derivatives can lead to 
internal inconsistencies that can cause a “locking” problem. As the length of 
the element goes to zero, the polynomial shape function approximation 
of the independent variable tries to match the approximation of its first 
derivative, which is a constant, when ‘locking” is present. This leads to the 
oscillatory behavior described above. But as more and more nodal points are 
added, i.e., more finite elements are added, oscillations of the independent 
variable will still persist if “locking” is present. These oscillations disappear 
for the STEP ID finite element as the mesh is refined and show that “locking” 
is not present for the one- dimensional cases studied here. 

To summarize, the depth of the boimdary layer mesh is determined 
by the larger of Ag// and ^(T). This is the starting point for determining 
the boundary layer depth. Numerical experiments are then used to expand 
the boundary layer until oscillations of |^1 disappear. Finally, additional 
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elements Eire inserted into the boundeiry layer imtil oscillations of Az also 
vanish. 

8.2.5 FIELD RECOVERY. 

The primary quantity of interest is again B. For this problem, there 
is no analyticcd solution, but the results can be checked to determine if 
they are physically correct. There are also two methods by which B can be 
determined from the finite element solution. Comparison of the results of 
these two methods determines if an internally consistent solution has been 
reached. 

The first method of determining the B$ field is the finite element 
approximation of Equation(5.3.2). The second method uses Equation(3.2.16) 
inserted into Equatition (7.1.5). The one dimensional form of this equation 
that gives the value for Bq at the outer node of each element e is 

S I rdrj (8.2.21) 

where the superscript letter in parentheses represents an element niimber. 
The integration over each element is performed by numerical quadrature 
using a two point Gauss rule. 

8.2.6 TEST PROBLEM. 

The test material used for this example was high purity aluminum. 
The materied constants a and for each element were evaluated at T equal to 
zero degrees Kelvin using the formulas of Section 4.4. The permeability fi of 
each element was set to fio as discussed earlier in this chapter. The reference 
state of V was set as described in Section 8.2.2. The mesh was discretized 
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as desribed in Section 8.2.4 with a regular mesh of Nbulk elements in the 
region 0 > r > rg — 200Ae//. Another regular mesh of Nbound elements was 
generated in the region — 200Ae// > r > Vc- N^uik and Abound denote 
the number of elements in the bulk and boimdary layers respectively. Nbuik 
and Nbound were 2 and 98 for this problem respectively. Because the free- 
space magnetic field element has been validated previously, all elements were 
within the conductor. The conductor radius Tc was 1.15 x 10 The value 
of lo wzis 5.0 Eimperes and the value of II was 0.0 amperes for the results 
presented here. The choice of these values ensures that an actual specified 
current loading for results presentation was attained. The element has been 
tested many times by loading from zero to full load and has worked extremely 
well. The only problem that was experienced was when the current loading 
appoached a magnitude that was large enough to move the solution close to 
the stationary point. In this region, Av and AC became increasingly smaller. 
To rectify this problem, the coding was modified to ensure that the step size 
at step n + 1 does not fall below an axbitr£iry value. If the step size became 
too small, .9 X In X Ii was added to J^, the reference state was reset, and 
the solution proccedure was restarted at step n. The only disadvantage to 
this scheme was determining the correct value of I at each step for output 
prirposes. This problem was easily circumvented by outputting the value for 
lo when it changed, and the step where the change took place. 

The main disadvantage with using the incremental solution methods 
for results presentation was that the solution process did not always stop at 
the desired full load value, but usually exceed it by some fraction of the step 
size In- This is a consequence of the solution procedure used, and is inherent 
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in arclength schemes. By setting II equal to zero and equal to J, this 
problem was bypassed here. 

For the results presented below, the solution procedure required 9 
iterations to converge, with the solution tolerance r being 4 x 10“^^. The 
2-norm for r did not include the value for the constrained degree of freedom 
of 1^|. The value of r there ranged between 10"^ to 10”’’ depending on how 
close the finite elements came to modeling a zero slope for \xj;\ at the conduc- 
tor /free space boundary. To more accurately match this slope requirement, 
all that was needed was a more refined mesh for this area. However the 
results obtained were judged accurate enough for our purposes. Finally, the 
estimated condition number of the system was 228. 

Figures 8.2 and 8.3 show the results obtained for the normalized val- 
ues of \ip\^ plotted over the whole conductor and the boimdary layer of the 
conductor respectively. If a London type superconductor had been modeled, 
an exact step function would have been expected. Because aluminum is an 
extreme type I superconductor, |0| should exhibit behavior that is almost 
“step”-like. Figure 8.2 illustrates that the finite element does model phys- 
ical behavior by returning values that closely match a step function. The 
boundary conditions are seen to match well in Figure 8.3 in that the slope of 
IV’Np is zero at the interior boundsuy and is very close to zero at the exterior 
boundary. 

Results for Az are shown in Figures 8.4 and 8.5. Figure 8.4 shows the 
behavior over the whole mesh and Figure 8.5 the behavior in the boundary 
layer. The physical behavior of Az should approximately be the opposite 
of |^|. Over the bulk of the conductor, Az should be zero, and where \tp\ 
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Figure 8.2: vs. r, values for complete mesh plotted. 
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Figure 8.3; l^p/jt/^ooP vs. r, values for boundary layer plotted. 
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Figure 8.4: vs. r, values for complete mesh plotted. 
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Figure 8.5: vs. r, values for boimdaxy layer plotted. 
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decreases, the magnitude of Az should increase as kinetic momentum is ex- 
changed for magnetic field momemtum. Figure 8.4 shows this expected phys- 
ical behavior. Figure 8.5 shows this behavior in more detail, and illustrates 
one difference in behavior between Az and |^|. The slope of Az is zero at the 
interior edge of the boimdary layer, but nonzero at the exterior edge. This is 
expected because the boundary conditions for Az and \ip\ are different at the 
exterior edge, the behavior of Az matching its expected physical behavior. 

Figures 8.6 and 8.7 display the results for jz over the entire conductor 
and the boundary layer respectively. The behavior of jz can best be described 
by making an analogy to a similar problem in fluid mechanics. The medium 
of the problem would be a large pool of water contained between two infimtely 
long straight walls. For convenience, the walls are aligned so that one is on 
our left side and the other on our right. To make the analogy correlate to the 
results presentation, the left wall would be the center of the superconductor, 
and the right wall would be the conductor/free space boundary. The bottom 
of the pool would be shaped so that the density of water molecules matches 
the density of the superconducting charge carriers. The walls and the bottom 
of the pool would present no resistance to water flow. Assuming laminar flow, 
a rapidly moving stream of water is injected into the pool along the right wall. 
For the EM problem, jz is analogous to the velocity of the water molecules, 
Uu, in the pool. Where the stream is injected, it is expected that a large, 
rapid change in Vy, would exist, which upon first examination would appear 
to be a Dirac delta frmction. 

This behavior is exactly matched by the velocities of the supercon- 
ducting electron pairs the finite element model, and is shown in Figure 8.6. 
At the conductor/free space boundary, a Dirac delta-like “spike ’ appears for 
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which is zero otherwise. A closer examination of the fluid velocities for 
the imaginary example wovJd reveal that would rise rapidly on both sides 
of the stream, but a more gradual rise in v^, would occur on the side of the 
stream facing the left wall as momentum is exchanged with the other water 
molecules there. Again mimics this behavior as show in Figure 8.7 and 
validates the ability of the STEPlD element to model the expected physics 
of a superconductor. 

This comparison of a fluid flow to a Ginzburg-Landau supercon- 
ductor is a particularly enlightening one because, if this superconductivity 
model is correct, it explains the source of the miniscule resistivity in super- 
conductors. The resistivity is a result of a momentum exchange produced by 
collisions of Cooper pairs, the supeconductor’s charge earners, as required 
by the residual equation (3.2.6). Because the collisions are relatively infre- 
quent, a “spike” in the current density appears in the boimdary layer, rather 
than a “smearing” of the current density to an approximate step function. 
The position of the spike is determined by the density of charge carriers, the 
cmxent density vector choosing the point where the fewest collisions can take 
place. The fact that the density of Cooper pairs is higher on the interior of 
the boundary layer than the exterior explains why changes more slowly 
towards the center of the conductor. The Cooper pairs in the current stream 
j'j are simply experiencing more collisions with stationary Cooper pairs be- 
cause the density of pairs is higher towards the center of the conductor. Its 
position also determines that there is an expulsion of the B fleld from the 
interior of the conductor (the Meissner effect) because there is no current 
there to generate a field in accordance with Maxwell’s equations. 
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Figure 8.8: B$ vs. r, values for complete mesh plotted. 
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Figures 8.8 and 8.9 show the Be field generated there. Figure 8.8 
shows Be over the whole conductor and Figure 8.9 shows Be in the boundary 
region. Because no analytical solution is available, the B field has been 
plotted using the two different methods cited previously and finite element 
values for |^| and A.^. As in Chapter VII, the Be field calculated using 
Equation (5.3.2) is plotted as a step function. Both sets of values match 
fairly well over most of the region, but show some divergence towards the 
maximum and minimum values of Be- No reason exists to prefer one set of 
values over the other, but using ^ to recover Be has the advantage of being 
able to directly compute Be at element nodes. 

Expected physical behavior is matched by both curves. The value of 
Be computed by using (8.2.21) also matches the necessary analytical value, 
derived from an integral form of Maxwell s equations, of figl A com- 
parison of these values with values obtained by using the London model of 
superconductivity does not allow any statement to be made about the accu- 
racy of the Ginzburg-Landau model because the former neglects the gradient 
of The important point is that the Ginzburg-Landau model must achieve 
a specific magnitude at r^., and this is verified by Equation (5.3.2). For the 
above reasons, the London values are not compared to the finite element 
values obtained here. 

8.3 FURTHER DISCUSSION OF RESULTS. 

A word of caution is necessary here with regard to the author’s phys- 
ical interpretation of results. Because no analytical solution is available for 
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comparison with the numerical results, these results and their physical in- 
terpretation must be treated with some suspicion pending experimental ver- 
ification. However, there is a good evidence to suggest that the results are 
valid. 

First, the numerical approach has been based upon the Ginzburg- 
Landau theory of superconductivity. This theory, while not being thoroughly 
validated experimentally for cases away firom the critical temperature, has 
been able to predict superconducting phenomena with a great deal of accu- 
racy (Ref. [21, pp. 104-191]). This provides a great deal of credibilty to the 
ability of the Ginzburg- Landau theory for the prediction of EM and quan- 
tiim phenomena within a superconductor. It is universally accepted eis an 
accurate model of the macroscopic quantum-mechanical and electomagnetic 
properties of a superconductor near its critical temperature 7^. 

Second, the results of this chapter and Chapter XI exhibit behavior 
that is in qualitative agreement with the physics of superconductors. These 
behaviors are the appearance of the Meissner effect and a current carried 
at the surface of a superconductor. The Meissner effect is the almost total 
expulsion of a magnetic field from the interior of a conductor, and this be- 
havior is shown in Figures 8.8 and 11.12. The cause of this effect is that the 
current is carried at the surface of the conductor (Ref. [31, p.335j). In order 
to satisfy Maxwell’s law Vx VxB = j, no current can be carried within the 
bulk of the conductor or a magnetic field will be present there. Again, the 
STEP ID finite element shows this behavior in Figure 8.6. 

Finally, there is some quantitative agreement between the STEPID 
finite element and a known physical value. The value of the Be field at the 
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conductor radius is known to be ijlqI / 27vrc* As mentioned at the end of the 
previous section, the finite element model achieves this value at 

8.4 SUMMARY. 

In this chapter, a broad range of topics necessary to the solution of 
the superconductivity problem by the finite element method are discussed. 
The topics include the four-potential formulation of superconductivity, ap- 
propriate boundary conditions, nonlinear solution techniques, scaling tech- 
niques, and appropriate mesh choices for fimte element models. The most 
important aspect of this research is the insight that is gained about super- 
conductivity. For the Ginzburg-Landau model, it is possible to think of 
the current that moves through a superconductor as a stream of charge 
carriers called Cooper pairs that moves through a “sea” of static Cooper 
pairs. This “sea” acts like an extremely low viscosity fluid, and the “stream” 
moves through the region of the “sea” where the density of the Cooper pairs 
is the smallest. This region represents the place where the least amount of 
energy is expended by the collisions of moving Cooper pairs with station- 
ary Cooper pairs. Unlike the London approximation, or linearized forms of 
the Ginzburg-Landau model, the physics of the system as described above 
are shown only by modeling the exact Ginzburg-Landau equations so that a 
complete description of jz can be obtained. The STEP ID model shows this 
behavior well, and from the limited search of literature that the author has 
performed, it is believed that this is the first model that shows the physics 
in such good detadl. 

Now that reasonable models for the normal and superconducting 
states of a conductor have been developed, the next step in the complete 
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modeling of a conductor is to add themiEil effects. This is the topic of dis- 
cussion of the next chapter. 




CHAPTER IX 


THE THERMAL PROBLEM 

It is first necessary to model the temperature distribution within a 
conductor before the effects a temperature field has on the EM fields and the 
quantum properties of a conductor can be determined. Appropriately, the 
first topic of discussion in this chapter is the modeling of the temperature 
field of the steady state heat conduction problem with convection cooling 
botmdaxy conditions. The one-dimensional case is the case of interest for 
this work's examples and is the only case discussed. In Chapter IV it is 
mentioned that there are no temperature gradients within a one-dimensional 
steady state superconductor. Because no gradients are present, the temper- 
ature distribution of a superconductor is known and the calculation of the 
temperature distribution by finite element methods is not necessary. There- 
fore, this chapter is concerned with the finite element modeling of a normal 
conductor. The temperature distribution within a conductor is a function 
of the current I and the thermal boundary conditions at Tc. In the current 
chapter, it is assumed that the current I is steady and does not change. 
Cases where the current load I changes axe discussed in the following chap- 
ter. For this chapter, the discussion is about the physics of a conductor as 
the thermal boimdary loads are varied. 

The discussion begins by first developing the finite element model 
for the temperature distribution of a one-dimensional conductor, and then 
determining the analytical solution of that problem. The zinalytical solution 
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to this problem is developed later because certain assumptions about the 
finite element model have a direct effect on the analytical solution. 


9.1 FINITE ELEMENT DISCRETIZATION. 


9.1.1 CONSTRUCTTNO THE LINTID FINITE ELEMENT. 

Using the two-node “line” finite element again provides the C° conti- 
nuity required by the variational functional of (4. 1.7) for T . Again individual 
elements and elemental properties are identified by the superscript (e). The 
two element end nodes are denoted by the subscripts i and j . The tempera- 
ture 7” is interpolated over each element as 

T = (9.1.1) 


where the row vector N contains isoparametric shape functions for the inter- 
polation of T. The elements of N £ire fimctions only of the spatial coordinate 
r. The column vector contains the nodal ^^ues of T, which are con- 
stants with respect to time. Substitution of these finite element assumptions 
into the variational fvmctional of Equation (4.1.7) gives 


7v(«) [2 Ot dr 

Variation of the above with respect to will produce 

^ JvM 1 




.(•> 


.(•> 


(9.1.2) 


(9.1.3) 
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Quantities k and are functions of the spatial coordinate r and not of the 
independent variable T, This assumption is made because the variational 
functional does not give the correct residuals for the heat conduction 
problem if the thermal conductivity and the electrical resistivity are allowed 
to be functions of T. This approximation can be corrected by using the 
nonlinear solution procedures of Chapter VIII. During the solution phase, k 
and are held constant at their values for step n, and after the solution 
vector V at step n + 1 is determined, k and are updated using the new 
temperature distribution of step n + 1. If the step size is small, v does not 
move too far from the true equilibritim path and the values of are close 
enough to the exact values that any error is negligible. This assumption is the 
reason for discussing the finite element model first instead of the emalytical 
solution. No analytical solution exists for the set of coupled EM-thermal 
equations where k and a; are functions of T. By maJdng the approximation 
that k and u are functions of the spatial coordinate r, an ansdytical solution 
can be developed. 

The spatial approximations used are the ones discussed in previous 
chapters, i.e., u is a step function or constant over cin element, and k is 
interpolated linearly across the element. In terms of our shape fimctions, k 
may be written as 

k = (9.1.4) 

where contains the nodal values of the thermal conductivity. The values 
for k^®^ are obtained by using Equation (4.2.1). Equation (4.2.1) is a function 
of T, and is evaluated at and to obtain the respective components 
of k^*\ k\^^ and For the evaluation of k^*^ and and 
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axe the two components of at step n. The elemental value for is 
determined by the formula 

u;('> = ^ (9.1.5) 


where Wo is the residual resistivity, and o;,(T) is the ideal resistivity. The 
ideal resistivity is calculated as described in Section 4.3. 

All but two of the boundary integrals of (9.1.3) can be ignored by 
noticing that the heat flux Qr between interelement boundaries in the radial 
direction must be continuous; this point is shown later by the analyt- 
icad solution. The two remaining boxmdary integrals occur at r equal to 
zero and At r equal to zero, it can be seen that the boundary integral 
there vanishes, and the only remsuning boimdary integral occurs at Tc- The 
remaining integral is a function of the boundary heat flux loading where 
Qr(C^) represents this load and is the thermal loading paramter. 

Because non-linear solution techniques are used, expressions for , 
f(®), and must be determined. The Euler equations of (9.1.3) 

give the following expressions for and f^®^ 


p(®) = 8m. {2Trr,Qr{C)} | { J | 

j.(e) _ f(e) _ p(e) 

where the Kronecker delta <5,j is defined as 

= 1 i = j 


(9.1.6) 


(9.1.7) 
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eind m = N^uirt + 1) with N^ire again representing the number of elements in 
the conductor. Talking the paxtials of r with respect to and produces 
the following expressions for and 


K(®) 




(9.1.8) 


where and u;(®) have been assiimed to be only functions of the nodal 
position r and therefore vanish. 


9.1.2 APPLYING BOUNDARY CONDITIONS. 

In the previous section, expressions for and q^®^ aire determined. 
They contaiin the boundary heat flux term C?r(C^)- The heat flux at for 
the one-dimensional steady state convection cooling problem is expressed by 
the Euler equation of (4.1.10) 

Qr = hconv{Too-nrc)) (9.1.9) 


Because the free stream temperature 7^ is known, it is natural to choose 
this \'ariable as the load to be varied. To make 7^ a variable load, it is 
split into two parts where 7^ is the initial loading temperatmre, and 7 l is the 
variable loading temperature. The free stream temperatmre is now expressed 
as 

= % + (9.110) 

The \'alue of the temperature at r equal to Vc can also be expressed as 


T(r,) = NjT; 






TTl — iV yntre "b 1 


(9.1.11) 
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where Nj is the second component of the shape function vector N, 
is the nodal value of the temperature at Tc, and the superscript letter in 
parentheses represents an element number and not an exponential. When the 
new expressions for Too ^(^c) substituted into (9.1.9), the expression 
for Qr(C^) becomes 

C?.(c’') = (t, + C’’7L - (9.1.12) 

and the expression for 5Qr(C^)/^C^ is 

= A„„„(Tx,-u,r>) (9.1.13) 


where is the nodal value of the incremental rate of change of T at Tc. 

Bv convention, any terms of a set of fimte element equations that 
include the incremental rates and the nodal displacements are usually moved 
to the left hand side of the system of equations. Doing this, eind making 
substitutions for Qr{(J) iii tli® previously derived vectors and gives 


0 0 
0 1 




f(<=) = J ^X^*^|+^me{27rr,} 

pi') = ^dvi') {u.I'IjW'n’’} + s„, (r„ + c’'Tl)} I ° I 


(9.1.14) 


Similarly, and are 
= / dV^‘^ ^ k 


I 

q^*^ = Sme {^“^^c^convTj} ^ 1 ^ 


dN'^dN 

dr dr 


d" ^me ■{27rrc) 


0 0 
0 1 


(9.1.15) 


The addition of the boundary term Qr makes this system of equations, when 
assembled, determinate and no nodal values need to be constrained. 
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9.2 ANALYICAL SOLUTION TO THE TEST PROBLEM. 

Equation (4.1.14) states that the temperature distribution for a lin- 
eeirly interpolated k and step function u and jz may be expressed as 


ujjz ( rAr kiAr^ 


Ak 


where Ci and C 2 are integration constants and 

, , kj — ki , Ak 

k = ki -\ — r — ki + — — r 


rAr 


Ar + Ak^ 


Tj- Ti 


Ar 


+ C2 

(9.2.1) 


(9.2.2) 


where ki and kj are the values of k at the inner and outer boundaries of 
integration respectively. To adapt this solution for each element, the inte- 
gration constants Ci and C 2 are first replaced by the constants and 

C^'^ven where the superscript (e) represents an individual element number 
again. For notational convenience, the following equalities are also defined: 
= k\^\ 6^*^ = AkfAr and . The function is also 

defined as 



(9.2.3) 


Using the new notation, the temperature T over each element is expressed 
as 

T(®)(r) = /i(®)(r) -h U(®)even (9.2.4) 


Using a little physics, it can be seen that the heat flux out of an 
element must equal the heat flux into adjacent elements because the system 
is conservative and energy must be conserved. This requires that the heat 
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flux Qr be C° continuous and gives the following series of equations that 
relate the heat flux transferred across adjacent element boimdaries 




C^^ldd ^ M 

.(«) 2 2 

j 


(9.2.5) 


This equation can be rearranged to give 

CHdd = ^^^+C('-^ldd- 




(9.2.6) 


2 ' ' — 2 

by using the relation = r\‘\ At r equal to zero, C^Hdd must also be 

for the system of equations to remain boimded. This gives for C^®)>dd 

c<‘U = ■^^-E‘=.‘>{¥('i''-'fO} 

= 0 


zero 


e > 1 
e= 1 


(9.2.7) 


The values for C^'lven may be derived in a similar manner by ensuring that 
T is C° continuous across element boimdaries. The result is 


C<‘>v.„ = Z- /t'-’C-S'*) - £' (9-2-8) 

P=« 

where 7^ is the temperature at the surface of the conductor and equals the 
nodal value Tj”^\ Tg is determined by using the discretized i-ersion of Equa- 
tion (4.1.21), which is 


Z = % + C^Tl + 



(9.2.9) 


nonzero 


The preceding analytical 
. As approaches 


solution is only \-alid for cases where 6^®^ is 
, the first two terms of (9.2.3) diverge. For 
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cases where 6^®^ is zero, the methodology used to derive (9-2. 1) in Chapter 
IV can be used to determine that becomes 

/(<<)(r) = + C<*Ud ln(r) (9.2.10) 

The methods of this section can then be used to show that the solutions for 
and C^‘\ven aJ'e identical to those of Equations (9.2.7), (9.2.8) and 

(9.2.9). 

The above solutions to the heat conduction problem were originally 
loaded into a Fortran subroutine and solutions to the excimple problem were 
computed. For the example problem and material values used, it was found 
that incorrect solutions were being determined. One source of error w^s that 
the magnitude of the first two terms of (9.2.3) were much greater than the 
last term. Finite precision numerics caused the last term to be virtually ig- 
nored when determining although this term should have 

made a noticeable contribution to both. All of the formulas were rearranged 
so that and were computed in a term by term manner, i.e., 

first all of the In(fc) terms were computed, then all of the ln(r/^:) terms, 
etc.. This improved the solution marginally, and the problem weis examined 
further. The largest source of error came from the finite precision mathe- 
matics again. The term 6^®^ was seen to be extremely small and caused the 
first two terms of (9.2.3) to diverge. Although the divergence of individual 
terms should cancel when siimmed during the computation of it was 

beyond the machine’s capability to resolve the minute differences between 
the large individual terms. False zero values or random values were being 
assigned for the difference by the machine. To correct this problem, when 
the absolute value of the percentage difference between and dropped 
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below 2.2 X 10~®, k was assumed constant over an element. The value of 
was then used for the elemental value of k, and Equation (9.2.10) was used 
to calculate and C^'iven- 

The numerical results became much better, and both of the above 
corrective procedures are implemented in existing coding. Results presented 
in this work as the analytical solution to thermal problems used the above 
methods to control numerical errors. 

9.3 NUMERICAL EXPERIMENTS. 

9.3.1 THE FINITE ELEMENT MODEL. 

The finite element model derived in the previous section has been 
applied to the test problem described later in this section. The LETID finite 
element is used for determining EM quantities and the LINT ID element is 
used to determine the temperatmre distribution. Both of these elements are 
treated els one- dimensional axisymmetric elements. The LETlD element is 
identical to the CUPLEID element except that is allowed to change 
during the solution process. The description of the nodal degrees of freedom 
and the variables associated with each degree of freedom for the CUPLEID 
element can be foimd at the beginning of Section 7.3.1. The permeability 
the resistivity and the current density are uniform over each 
element. 

For the LINTID element, the “line” type element has only two end 
nodes which axe defined by their axial positions and They each 

have one degree of freedom corresponding to the temperature T which re- 
sults in a total of two degrees of freedom per element. These nodal values are 
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determined by interpolation with standard linear shape functions that pro- 
vide the C° continuity required by the veiriationcil formulation. The thermal 
conductivity of each element is also calculated by interpolation with stan- 
dard linear shape function where the nodal values of k are determined by the 
nodal temperatures at step n and use of Equation (4.2.1). The resistivity of 
each element, for both the LETID and the LINTID elements, is calculated 
by using Equation (4.3.1) and the nodzJ temperature values at step n to 
determine w at each node, and then taldng the mean of the two a;’s. The 
value of for the LINTID elements is determined by use of the LETlD 
finite element. 

9.3.2 APPLYING BOUNDARY CONDITIONS. 

As shown earlier, no nodes are constrained for the thermal part of 
this problem. The thermal flux terms that contain the Kronecker delta are 
directly injected at their appropriate positions when assembling f^®\ 
r(e)^ and to accoimt for boundary conditions. The electromagnetic 
boundary conditions are set as described in Section 7.3.2. 

9.3.3 ASSEMBLY AND SOLUTION. 

Both tangent stiffness matrices and are assembled in an el- 
ement by element fashion following standard finite element techniques. 
and K are used here to represent the master electromagnetic and master 
thermal tangent stiffness matrices respectively. The superscripts EM and T 
are also used in the sequel to distinguish between assembled electromagnetic 
and thermal vector quantities (e.^., is the electromagnetic solution vec- 
tor). is stored in a symmetric skyline form and is stored as three 
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vectors because it is a tridiagonal matrix. Systems of electromagnetic equa- 
tions are modified for boundary conditions as described in Section 7.3.2 and 
are processed by a standard symmetric skyline solver. After the solution pro- 
cedure that is described later has been used, provides the desired nodal 
values of A 2 for field recovery and the elemental values of jz for calculation 
of the electromagnetic heating loads of the heat conduction problem. 

Svstems of thermal equations are modified for boundary conditions, 
as discussed in the preceding section, and are then processed by a standard 
tridiagonal solver. The solution procedure then returns which contains 
nodal values of T . 

Solution Technique 

Because the values for k and are actually functions of the tem- 
perature T and not the spatial coordinate t , the LINTlD finite element is 
nonlinear, and the nonlinear solution techniques of Section 8.2.2 are used to 
solve problems. These techniques work well for thermal problems if the rate 
of change of temperature across an element is not too large. 

The solution procedure is started by choosing a reference state for T . 
The reference state chosen for the examples of this work is set by initializing 
to Tq and to zero. For cases where is not sufficiently close to 

the equilibrium path, the reference state may be chosen by use of Equation 
(4.1.22). For the latter case, the thermal conductivity and electrical resistiv- 
ity are evaluated as constants over the whole domain of the conductor and 
evaluated at some mean representative temperature such as %. The ciurent 
density becomes a constant with these asstimptions and is equal to 
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To advance the solution to the next step, the set of incremental rate 
equations 

K„w„ = q„ (9.3.1) 


is solved. This system can be written as 


K 


EM 

n 

0 


0 


K 


n J 


W 

W 


EM ■) ( „EM 


where the subscript n represents the current step nmnber. The solution 
vector v„ is equal to }^. The rest of the solution procedure is 

identical to the procedure outlined in Section 8.2.2. except that K, r, f and 
p at step n + 1 are calculated by using the values of T and from the 
previous step n. The solution must be calculated in this manner because 
the variational formulation assumes that k, w and axe all functions of the 
spatial coordinates tind not the independent variable T. This device holds 
the material properties and constant for K, r, f and p so that a new 
temperature distribution at step n + 1 cein be determined. This means that 
the solution procedure is not solving the correct set of equations for the true 
equilibrium path of the heat conduction problem. 

Unfortunately, this is the current level of axivancement of heat con- 
duction analysis. To solve this problem, a relatively small step size is 
tahen. This “fix” allows the computed solution to remain close to the cor- 
rect equilibrium path. This problem can also be corrected by holding 
constant eind taldng several steps. The result is that the quantities K, r, f 
and p are all updated until the correct equilibrium path is reached. For the 
examples presented here, the last “fix” was not required as taking a small 
step size /„ brought the solution sufficiently close to the true equilibrium 
path. 
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9 . 3.4 5^C,4LING TECHNIQUES. 

The tridiagonal solver that was implemented does not include any 
method for the estimation of the conditioning of a system of equations. Be- 
cause no estimate was available, no attempt was made to scale any of the 
thermal systems of equations that were processed. The skyline solver did 
include the capability to estimate the condition number of K . Choosing 
realistic values for the material properties and made so 

highlv ill-conditioned that the solution method failed, and showed that the 
values used in the example problems of Chapter VII were rather simplistic. 

Some of this ill- conditioning is alleviated by employing a scaling 
scheme similar to the first scaling technique of Section 8.2.3. This scaling 
is performed at an elemental level before is assembled. Using T, M, 

T and Q to represent units of length, mass, time and charge respectively, 
the units of Az, jz, and Ac axe seen to be ML/{TQ), Qf(TL^)j QIL and 
ML?' 1{TQ) respectively. Many scaling schemes were tried, but the one that 
reduced the condition number the most gave the scaled displacements of Az , 
jc, and Ke dimensions of M^I'^LfT and Ac dimensions of IT. The 

elemental scaling matrix for each element is 


c(e) 

“Jll 

0 

0 

0 

0 

0 

0 


c(*) 

^22 

0 

0 

0 

0 

0 


c(*) 

‘->33 

0 

0 

0 

0 



c(*) 

“^44 

0 

0 

0 






0 

0 



symm. 



c(e) 

•^66 

0 







c(c) 
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o(«) 

“^11 

c(«) 

— O44 

= 2/^ 

c(«) . 
<^33 
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c(«) _ 
^22 ~~ 
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c(e) 

*^55 

= 0 e 
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where is the elemental length The inclusion of the elemental 

length in the seeding parameter for jz does not cause any difficulties because 
the only diagonal terms affected by 5^2^ of correspond to the jz degrees 

of freedom which do not couple through the diagonal terms. 

To illustrate this point, consider the diagonal terms corresponding to 
the degrees of freedom for Az of a two element system that contains a toted 
of three element end nodes. K A^z^ is scaled for each element by a[^^ at 
the shared center node will be scaled on an elemental level by for the first 
element and for the second element. Note that the bracketed superscript 
numbers represent the element number and not an exponent. This causes 
no problem in determining scaling factors if equals It is simply 

or But if the lengths differ, problems will occur because Az at the 
shared node will be scaled into two different dimensions! Trying to assemble 
the scaled elemental diagonal terms of Az with this scaling would result in 
an error because each scaled variable represents a single scaled independent 
variable. 

The scaling matrix of (9.3.3) avoids these difficulties due to a careful 
choice of its elements. ensures no coupling of adjacent values of and 
also ensures that there are no zero diagonal terms of the assembled scaling 
matrix for degrees of freedom of that are constrained to zero. also 
ensures that no zero diagonal term appears for the extra “empty” degree of 
freedom of that is discussed in Section 7.3. 
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The scaled elemental electromagnetic tangent stiffness matrix for 

1 ^ e ^ ^ xvirc 


ic"' = ^ 


/(e) 


0 


(e) 

-r) 

* m 


symm. 


A^23 = 


/(e) 




2 w(*) 
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A^23 


0 
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0 

0 




(e) 




Tm 

/eV 

A^26 


o . J 


0 

K27 

0 

0 

0 

0 

0 


(9.3.4) 


fioCo _(e) 

where = (l/2)(r-*^ + Tj*^) is the meein radius of the element and Cg 
is the speed of light in vacuum. As can be seen by inspection, the closer 
(/^®V2)(/^oCo/w^®^) is to one, the better the conditioning of the matrix. If 
this occurs, it can eiIso be seen that all of the terms become proportional 
to approximately Using the test material of high punty almninum 

does not allow to come as close to one as is desired, and 

other means are used to reduce the conditioning of the scaled electromagnetic 
system of equations. It was found that the choice of mesh discretiztion 
greatly affected the conditioning of the system, and this is the next topic of 


discussion. 
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9.3.5 ME8H GENERATION. 

The thermal conduction problem is similar to the superconductivity 
problem because a finer mesh discretization is required neair Tc to accurately 
determine the nodal values of the independent vciriables. The independent 
variable for the thermal conduction problem is T cind the electromagnetic 
quantity u in the temperatxire range of interest is proportional to . This 
behavior suggests that electromagnetic quantities vary either much more 
quickly or more slowly than T depending on whether the system is above 
or below 1° Kelvin. This behavior also suggests the use of separate meshes 
for the electromagnetic and thermal equations. Separate meshes were not 
used for the examples of this work because they require a transfer of data 
between the thermal mesh and the electromagnetic mesh and are subject to 
extrapolation errors. Separate meshes also require more computational effort 
and memory storage. For the above reasons, it was decided to use a single 
mesh for both the linear electromagnetic and thermal system of equations. 
As mentioned above, there is a finer grading of the thermal mesh at r^. This 
grading was determined to be a source of ill-conditioning for the assembled 
EM equations. 

Originally, a small region near Tc was discretized with regular finite 
elements and the remainder of the conductor volume was descretized with 
larger regular finite elements. It was seen that at the node where the two 
meshes joined, off-diagonal terms were generated that were substantially 
Itirger than diagonal terms. Other off-diagonal terms that were substantially 
smaller than diagonal terms were also generated. To cure this problem, the 
seeding scheme of the previous section is used to make off-diagonal terms of 
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th.e same order of magnitude as the diagonal terms on an elemental level. A 
new mesh discretization scheme is also used to eliminate the substantially 
larger and smaller olf-diagonal terms. The third scaling technique of Chapter 
VIII does this for a system of positive definite equations, but this technique 
was tried and did not work here because the system is positive semi-definite 

due to the Lagrangian multiplier Ac. 

By observing that all of the terms in the elemental EM matrices are 
approximately proportioned to basis for the mesh discretization 

can be determined. By minimizing the rate of change of bad 

off-diagonal terms can either be eliminated or changed to more closely ap- 
proximate the magnitude of the diagonal terms. 

For convenience, a regular mesh of Nfine elements is used in the 
region near Tc. The remaining region uses a special mesh that contains 
Ncoarse elements. This choice determines d(rlt ? to be 1.0 in the 
“fine” region. It also specifies that equals (rncpi/l^*^)"k -b a-t e equal 

to Ncoarse + 1 where r„cpi is the value of r at the node where the coarse 
mesh ends. This choice also determines that the length of each element in 
the “fine” region be equal to (rc — r„cpi)IN fi„e- An additional boimdaiy 
condition for the axisymmetric problem is that rm ^ is always equal to .5 
at e equed to 1. 

To satisfy the three boundary conditions and still have a varible 
input for the mesh discretization, a cubic curve fit for r\^ is required. 
Using these requirements gives the following equation for 

(e) ,, 

^ = I + A^e'^ + A^e -h A.) (9.3.5) 

Ai = (2Ci — BiNcoarse {Ncoarse ~ 1) “ ^C2{N coarse + 1) + ^coarse ~ 1) 


Kill I iiiiiii mill II I I II nil II iiiiii i 1 1 ii .iiiiiiitii III! I in II mil 
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A2 = Bi-^ B 2 A 1 Az = 2 C 2 - lAi - 3^2 - 1 


Ai = 1 — (^1 + A2 + 


2C3 — 2C2 + 1 


Bi = 

Bz = (N^oarse + 1 )" + B2{N, 


coarse 1 
\3 


B 2 




— 3AT^ 

coarse 


2Nc 


- 1 


+ 1)^ — (7 + ^B2^{Ncoarse + 1) + 2^2 + 6 


Cl = + .5; e = Ncoarse + 1 C2 = e = 2 

where the subscripted ^’s, 5’s and C’s represent constants. The constant 
C2 is an input value amd represents the value of fl^‘^ for the second el- 
ement. The value of for the second element is chosen as in input 

variable because at the first element, this value is determined by the prob- 
lem geometry as always being 0.5. It is also easier to determine the size of 
the second element with this formula. If C 2 is larger (smaller) than 1.5 the 
second element is smaller (larger) than the first, and if C 2 equals 1.5, it has 
the same size. The value for r where the two meshes meet (r„cpi) is also 
an input value. For the numerical experiments presented here, it was found 
that Tncpi = .75rc and C2 = 1.5 produced the best results. 

To determine the values of r at each node in the cozirse mesh region, 
the following additional formula is used 


(e) _ (0 ( . 


C — (i\r coarse “I" l)? coarse? • - - , 2, 1 (9.3.6) 
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9.3.6 FIELD RECOVERY. 

For the thermal part of the problem, the T field is recovered by- 
using the values of v^. The electromagnetic fields are determined by using 
the simple finite element approximation of Chapter V, Equation (5.3.6). This 
value was plotted as a step function due to its <7 ^ continuity. 

9.3.7 TEST PROBLEM. 

In this example, high purity aluminum was again used for the test 
material. The values for k and were calculated as discussed earlier in 
this chapter. The permeability and permittivity for all elements were set to 
fio and as explained at the beginning of Chapter IV. The reference state 
for V was set as described in Section 9.3.3. The geometry was that of a 
one-dimensional axisymmetric wire as shown in Figure 2.1 with a radius Vc 
transporting a total current I equal to 5 amperes in the positive z direction. 
The mesh was discretized as described in Section 9.3.5 with Tc equal to 
1.15 X 10“^, Ncoarse equal to 50 elements and Nfine equal to 30 elements. 
This gave a total of N^^irt equal to 80 elements. Because the element for 
the free space magnetic fields had been validated before, no elements were 
generated external to the condnctpr. The solution tolerance t was 9.0 x 10 
and required about 2 iterations per step to converge. The estimated condition 
number for ranged from 10816 to 65888. The initial temperature 'was 

chosen as I® Kelvin and the loading temperature w^ 1" Kelvin. The step 
size /„ was chosen as .025 and 40 steps were necessary to move the solution 
from the starting temperature of 1“ Kelvin to 2° Kelvin. 

The results for the analytical solution and the finite element solution 
for T were so close as to be indistinguishable on a plot. Consequently, Figure 
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9.1 shows the results for the finite element solution for the temperature 
distribution and Figure 9.2 shows the percentage error of the finite element 
solution from the analytical solution as functions of the radial distance at 
the final step (To© = 2° Kelvin). A comparison of the Euler equations of the 
EM problem (Equation (2.3.11)) and the thermal problem (Equation (4.1.8)) 
shows that they are of an identical form. If fc is nearly constant and 
approximates a constant, then the behavior of A* Bind T within the example 
wire should be the same. A comparison of Figure 9.1 with Figures 5.2 and 

7.2 show that this is indeed the case. Figure 9.2 shows the deviation of the 
computed solution from the cinalytical solution as a percentage error, and it 
ranges from a maximum at r = 0 of 2.860 x 10~^ to zero at r = r©. 

The primary variables of interest for this research were the B fields, 
and they are the only EM results that are shown here. Figure 9.3 and 9.4 
show the results for the Be field at the final step. Figure 9.3 shows the Be 
field over the whole conductor and Figure 9.4 shows the B$ field for the 
volume of the mesh with the finer discretization. The percent error ranged 
from 33.1 percent at r = 0 to 5.89 x 10“^ at r = r©. By observing where 
the Einalytical solution intersects the “steps” in Figures 9.3 and 9.4, a rough 
estimate of the accuracy of the finite element solution can be made. The finite 
element solution is exact when the analytical solution intersects the center 
of the “step tops” . It can be seen that the analyticzd solution intersects the 
majority of the “steps” at their center points. It can also be seen that the 
error quickly diminishes as the distance from the center of the conductor 
increases. 

As mentioned earlier, the solution procedure is not exact. Also men- 
tioned was that the exact solution can be computed by setting % to the full 
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Figtire 9.1: T vs. r, values for the finite element solution plotted. 
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load, 7 l to zero and using the corrective Newton-Raphson method to iterate 
onto the equilibrium path. 

The above technique was tried for this numerical example and it was 
fo un d that the finite element results presented here differed from the exact 
solution by 4.65 x 10~^ and 3.43 x 10”® percent for the temperature at r equed 
to zero and Tc respectively. The B$ field differed from the exact solution by 
33.3 and 5.60 x 10“'* percent at r equal to zero and Tc respectively. 

A brief word must be said here about the analytical solution of Sec- 
tion 9.2. In general, the analytical solution does not match the exact solution 
because it uses the nodal values of T from the previous step to compute fc, 
and The analytical solution only becomes the exact solution when 
is held constant and the solution is allowed to iterate onto the equilibrium 
path. This brings about the rhetorical question, why bother computing the 
analytical solution? 

The zmalytical solution is computed because it gives some measure 
of how close the finite element solution is to an exact solution. It will always 
give the correct form of the solution, but not the correct magmtude if the 
solution lies close to the true equilibrium path. This is the first step in 
assessing the accuracy of the solution vector produced by the non-linear path- 
following techniques used in this work. The second step is to hold the loading 
parameter constant and then use the corrective Newton-Raphson technique 
to iterate to the exact solution. The exact solution is then compared to the 
original incremental solution. If the difference between the two solutions is 
too large, then the solution procedure has moved too far from the correct 
equilibrium path md a smaller step size needs to be chosen. The solution 
process is attempted again with the new step size and reset to zero. 
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Figure 9.3: Be vs. r, values for complete mesh plotted. 
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Figure 9.4: Be vs. r, values for finer discretization plotted. 
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In general, it is not possible to hold the loading parameter constant 
at its final %-alue and iterate onto the exact solution without first following the 
incremental path to that point. There may be bifurcation points or other 
stationary points on the equilibrium path that will not eillow the solution 
method to converge to the correct solution by this simple iterative process. 
The problem of the thermally loaded conductor presented here does allow 
the above method to converge because it can be determined where the only 
critical point for this problem lies, that point being the superconducting 
phase transition point. 

The plots for T and presented here show that with an appro- 
priate mesh choice, a reasonable step size and a slowly varying temperature 
distribution that the solution technique presented here is adequate for the 
author’s purposes and little accuracy is lost with this solution procedure. 

9.4 SUMMARY. 

In this chapter, it is shown how thermal fields may be modeled with 
the LINTID finite element. The CUPLEID finite element is also adapted 
to the nonlinear solution techniques of Chapter VIII to become the LET ID 
element. This demonstrates the usefulness of the four-potential method and 
the solution techniques for modeling the coupling that occurs between ther- 
mal and EM fields. The four-potential theory is also validated for computing 
the desired EM quantity, the B fields and the effects of temperature on these 
fields by the results presented in this chapter. The use of real values for w, 
and e make the problem more difficult to solve, but by using a different mesh 
discretization and scaling techniques, good solutions for the B field can still 
be readized. 
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For this chapter, the case where only thermal loads are allowed to 
vary is solved. In the next chapter, the case where the EM load I is allowed 


to vary is examined. 


CHAPTER X 


COUPLED THERMAL-EM PROBLEM IN NORMAL CONDUCTOR 

In the previous chapter, the consequences of varying the ther- 
mal loading on a conductor eire discussed. In this chapter, the thermal- 
electromagnetic coupling in a normal conductor loaded by varying the cur- 
rent I is discussed. Most of the necessary ground work to consider this 
problem is developed in the previous chapters. Analytical solutions to both 
problems axe discussed in previous chapters and axe not presented here. The 
solution, mesh discretization and scaling techniques of the previous chapter 
are also implemented for the varying current load problem. The only paxts of 
this problem that change axe paxts of the LINTID finite element that depend 
explicitly upon j, which is a function of the current load I, and the paxts 
of the LETID finite element that are dependent upon I. The first topic of 
discussion is the modification of the LETID finite element to include cases 
where I is allowed to vary. 

10.1 FINITE ELEMENT DISCRETIZATION. 

10.1.1 MODIFICATIONS TO THE LETID FINITE ELEMENT- 

The first step in adapting the LETID finite element for a varying 
current load is to split I into an initial current lo and a loading current II 
as wais done for the superconductor. Thus I — !„ + where is 

the electromagnetic loading parameter. By using the Kronecker delta with 
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the boundary terms, the new equality for J, emd Equation (7.2.2), 
and are expressed as 

f(e) _ (10.1.1) 

*0 0 0 0 0 0 
0 0 0 0 0 
0 0 0 0 0 

= 2ttH 0 0 0 0 

^me 0 0 

symm. 0 0 

0 . 

( 10 . 1 . 2 ) 

f S^eH 
0 
0 

p(') = (/^ + < 0 

0 

V. ^me 

where m again equals iV^.ve + 1 and H represents the element height. The 
matrix is used here to add the boundary terms for k to ^ discussed 

in Section 7.2.2. also removes the rank deficiency generated by the 

“empty” degree of freedom discussed in Section 7.3.2. The correct value of 
is given by Equations (7.2.3), (7.2.4), (7.2.5) and (7.2.6). Taking the 
partials of and p^®^ with respect to the independent variables contained 
in gives the tangent stiffness matrix and loading vector, respectively. 
Therefore the tangent stiffness matrix K^®\ the loading vector q^®^ and the 
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incremental rate vector w 


(e) 


axe: 


K 




f S,.H -i 
0 
0 
0 
0 

-Sm.H 


> ; w 


(«) _ 



5C 


dC 


£m 


(10.1.4) 


It can be seen that is exactly the same as the stiffiaess matrix, modified 
for boundary conditions and the extra “empty” degree of freedom of Chapter 
VII thereby justifying the statement at the end of Section 7.1.1 that the two 
are equivalent. A word of caution is necessary here because the expressions 
for and are only valid for cases where the first degree of freedom for 
Az is constrEiined. For the excimples presented in the sequel, this is always 
the case and no further information is needed to solve the EM system of 
equations except that the “empty” degree of freedom for the fifth element of 
must also be constrained to zero as explained in Section 7.3.2. 
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10.1.2 MODIFICATTONS TO THE LINTID FINITE ELEMENT. 

The expressions derived in Chapter IX for K^, and for the 
LINTID finite element do not depend upon the current load I explicity 
and may be used without modification to solve the current loading problem. 
Because for the LINTID element does contain \ it must be varied 
with respect to the loading parameter to determine what happens to 
the thermal system of equations as / is varied. The result for is 


q 


(e) _ 



(10.1.5) 


The assembled system of coupled thermal and EM equations with this form 
for q can be expressed as 


K 


EM 


0 




where q^ is now a function of . This form of the equations is undesirable 
because terms of the incremental rate vector appear on both sides of the 
equality sign. Moving q^ to the left hand side produces 


kEM 

Y^emt 



(10.1.7) 


where contains the elements of -q^ at the appropriate positions. 

This form of the equations is also imdesirable for two reasons. First, 
the extra matrix must be assembled which requires more computa- 

tional effort and memory storage. It also ruins the sparsity and the symmetry 
of the original system of equations and a tridiagonal solver can no longer be 
used to process the thermal equations. Second, is also a function of 

the independent variable This affects adversely the conditioning of the 
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coupled system and typiccilly requires more iterations to converge upon a 
solution than a set of linear equations, which only requires a single iteration 
to converge. 

If the form of Equation (10.1.6) is preserved by making a function 
of , the original sparsity of the system is ret 2 iined and computational 
effort ajid memory storage are reduced. To accomplish this goal, Equations 


(7.1.2) and (7.1.3) zire used. Insertion of (7.1.3) into (7.1.2) gives 

ire A ,(P) 




i^uitre /* ,(P) 

p=i •'M 


where the superscripts (p) and (e) again represent element numbers and 
Njpire is the total number of elements within the conductor. Rearranging 
Equation (10.1.8) and taking the partial of with respect to gives 


1 

q^EM ^^(e) 


[ /r« ] 


(10.1.9) 


The bracketed term of the above expression is evaluated when assembling 
and requires less computational effort and memory storage than the 
scheme presented in Equation (10.1.7). The thermal loading vector is 
still a function of the EM solution vector but does not cause diflBculties 
in the example problems. The value of at step n is used to compute w 

at step n -r 1 in the same manner that is used at step n to compute k 

and u: for and at step n + 1. 

The coupled system of thermal and EM equations is conservative 
and the \'ariation of the discretized functionzJs that describe this system 
should produce a symmetric system of equations. The fact that (10.1.7) is 
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not symmetric is caused by the use of approximations that render them in- 
complete. In Chapter IX, it is assumed that j is only a function of I but not 
of the temperature T, an approximation responsible for the unsymmetry Of 
(10.1.7). This point is mentioned here to emphasize the importance of using 
the nonlinear solution procedure of Chapter VIII and equilibrium path fol- 
lowing procedures to generate correct solutions. For small incremental steps, 
the methods described in this and the last chapter work well for determining 
solutions for the coupled system of thermal and EM equations although the 
complete set of equations is not solved. Because the solution never moves 
too far from the equilibrium path, the missing terms have little effect on 
determining a correct solution. 

10.2 NT.'MERICAL EXPERIMENTS. 

10.2.1 THE FINITE ELEMENT MODEL. 

The finite element model described in the previous section has been 
applied to the infini te axisymmetric normal conductor of the previous chap- 
ter. The modifed LETlD and LINTlD finite elements are used to determine 
EM quantities and the temperature distribution, respectively. Both of these 
elements are treated as one- dimensional axisymmetric two node “line” type 
elements. The description of the nodal degrees of freedom and the vari- 
ables associated with each degree of freedom for the LETlD element are the 
S 2 une as the CUPLElD element and are found at the beginning of Section 
7.3.1. The permeability fx, the resistivity w and the current density jz are 
all eissumed to be uniform over each element. 
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For the modified LINTID fimte element, the description of the nodal 
degrees of freedom and the variables associated with each degree of freedom 
are found in Section 9.3.1. The material constants k and u; are also deter- 
mined by the same methods discussed in that section. The value used for 
to determine at step n -H 1 comes from the EM solution vector at 
step n and for is determined by use of Equation (10.1.9). 

The boundary conditions for this system zire the SEime as described 
in Section 9.3.2 and are set in the same manner. 

10.2.2 ASSEMBLY AND SOLUTION. 

The assembly and solution of the coupled EM-thermal normal con- 
ducuctor with a variable current loading is identical to a conductor with 
variable thermal loading except that Equation (10.1.5) is used to determine 
q^. The scaling techniques implemented for for the thermally loaded 

conductor are also used on The mesh generation and field recovery 

techniques used here are also identical to the techniques of Sections 9.3.5 and 
9.3.6 respectively. 

10.2.3 TEST PROBLEM. 

A one-dimensional axisymmetric wire made of high purity aluminum 
was used as test example. The geometry is shown in Figure 2.1. The radius of 
the wire Tc is 1.15 x 10““*. The wire transports a total current of -|- 
in the positive z direction where lo is zero amperes and Ii is 5 amperes. 
The free stream temperature of the system 7^ is held constant at 2° Kelvin 
by setting the initial temperature %, equal to 2" Kelvin md the loading 
temperature 7 l equal to zero. The convection heat transfer coefficient hconv, 
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the thermal conductivity k and the resistivity are calculated as described 
previouslv in Chapters IV and IX. The permeability and the permittivity 
for each element were set to fig and respectively, which are the free 
space or ■\'acuum vtdues, 2 is discussed at the beginning of Chapter IV. 

The mesh was discretized as described in Section 9.3.5 with Ncoarse 
equal to 50 elements and Nfine equal to 30 elements to give the total of 
elements within the conductor, Nwin equal to 80 elements. Because the 
element for the free space magnetic field had been validated before, none 
were used outside of the conductor. The incremental step size /„ chosen was 
.1 and 21 steps were taken to give a final value for of 0-97. The 
solution tolerance r was 1.0 x 10“^ and the solution procedure averaged 
2.381 iterations per step to converge. The condition number of was 

estimated to range from a low of 14826 to a high of 70969. 

The results of the analytical solution and the fimte element solution 
for T axe indistinguisable on a plot. Figure 10.1 shows the results of the 
finite element solution and Figure 10.2 shows the percentage de'V'iation of 
the finite element solution from the analytical solution for the final value of 
^EM Tiie maximum error from the zmedytical solution occured at r equal 
to zero and was 2.67 x 10~®. 

To converge upon the exact solution another four incremental steps 
using 14 iterations per step was reqtiired. For these four steps was held 
constant. The results of this final process showed that the finite element 
solution presented in Figme 9.1 differed by 9.36 x lO^^ percent at r equal to 
zero and_5.27 x 10“® percent at r equal to from the exact solution. For 
all increments, the initial or reference state was the finite element solution 
obtained after following the incremental path to equal to % 0.97. The 
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state at w 0.97 was actually quite close to the equilibrium path. The 
high iteration number required to move the solution back onto the eqioilib- 
rium path illustrates the difficulty associated with finding an exact solution 
by simply using the corrective Newton-Raphson process. Sometimes, direct 
iteration is not possible or more expensive computationally than just using 
the incremental “path following” solution method presented here. 

Figmes 10.3 and 10.4 show the finite element and analytical solutions 
for Bg. Figme 10.3 shows the two solutions plotted over the whole domain 
of the conductor, and Figure 10.4 shows the solutions over the domain of 
the finely graded mesh. At r equal to zero, Bg differs from the analytical 
2 ind exact solutions by 33.3 percent. At r equal to r^, Bg has errors of 
5.91 X lO”"* and 5.85 x lO”'* perecent when compared to the analytical and 
exact solutions respectively. 

Although not mentioned until this point, the 33.3 percent error at 
r equal to zero appears to be large. Subtracting out the temperature at r 
equal Tc from eJl of the finite element and analytical values for T and then 
recomputing the percent error for T will give the same error at r equal to 
zero of 33.3 percent. This observation has led the author to believe that the 
error is the same for all systems of equations of the form V • (aV&) = f{b) 
when this system is modeled with finite elements. Here a represents some 
material constant, b the independent variable, and /(fe) some function of the 
independent variable b. If the above observation is correct, it should always 
be possible to correct any error when modeling a system of this form. 

A quick look at Figures 5.3, 5.4, 7.3, 7.6, 9.3 and 10.3 shows that the 
correction is probably unnecessary. The divergence of the computed solution 
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from the exact solution is so small as to be almost unnoticeable and from a 
practical engineering standpoint, the error is not noticeable. 

The reason that the 33.3 percent error appezirs to be large is because 
the computed solution is compared to the exact solution on a node to node 
bzisis. This is formally expressed as 

% error = X 100 (10.2.1) 

where Bfe and Bex are the finite element and exact values respectively for 
the Bg field. A more realistic error estimator for engineering purposes is 


0 / BpEir) - BEx{r) _ 

% error = — r-r :: — r X lUU 


( 10 . 2 . 2 ) 


Bex(^) + Bex{^c) 

This type of error estimate has been used to compute the error for T in 
this and previous chapters. When computing T on a node by node basis, 
the boimdary loading has already been factored into the estimator. The 
conclusion of this brief digression is that an error estimator by itself is not 
always a good indicator of the accuracy of a finite element model. Graphics, 
a relationship of numericcd answers to the actual physics of a modeled prob- 
lem and good engineering common sense should all be used with an error 
estimator to judge the validity and usefulness of each finite element model. 


10.3 SUMMARY. 

In this chapter, a form of the LINTID finite element is derived for the 
c£ise where I is varied instead of T. The LETID finite element is modified, 
and using a nonlinear solution technique it is possible to compute some good 
values for the thermal and magnetic fields of a one-dimensional axisymmetric 
conductor. 




Figure 10.4: B$ vs. r, values for finer discretization plotted. 
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This represents the next to the last phase of this work s analysis of 
the coupled quantum phase- EM-thermal problem for superconductors. It 
is now possible to generate models of the superconducting material in its 
normal and superconducting phases. It is also possible to determine the 
effect of ■^'arying either the thermal T or current I loading of that material. 

The next chapter is concerned with the use of the final versions of 
the LINTID and LETID finite elements with the STEPID finite element 
in a single computer program. This program is used to determine the cor- 
rect state of a superconducting material and the values of the thermal and 
magnetic fields. 




CHAPTER XI 


THE COMPLETE COUPLED PROBLEM. 

This chapter is concerned with the correct application of the LINTID, 
CUPLEID and the STEPID finite elements to the the specific problem of 
determining the electromagnetic and thermal fields within a superconducting 
material. The main limitation on this model is that it is only one-dimensional 
and caimot realistically model the state where B and 'T reach their critical 
^•alues because, as noted previously, the solution at the transition state is a 
multi-dimensional problem. At this bifurcation point, a mixed normal and 
superconducting state appears that needs to at least include variations of 
li' in the z direction to obtain an accurate model of the physics that occur 
wdthin a conductor [21, pp.99-103; pp.127-191]. This point also marks where 
the computed solution must change equilibrium paths to accurately model 
the physics of the electron transport within a conductor. 

The methodology required to model a superconductor when the ther- 
mal loading is varied is first discussed in some detail. Then the question of 
how to determine if a computed solution lies above or below the bifurcation 
point is discussed. Means to determine whether or not such a solution is 
correct are also presented. These means form the basis of an equilibrium 
path changing criteria. Finally, results that show the model changing state 
as I and are varied are presented to validate the path changing criteria. 
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11.1 A SUPERCONDUCTOR WITH A VARYING T LOAD. 

To use the incremental approach to solve the superconductor with 
a varying thermal load, and are varied with respect to and 

The thermal quantities and are also varied with respect to 
Doing so produces a system of incremental equations that can be partitioned 


in the following manner: 


K 


EM 


K 


EMT-\ 


0 


K‘ 




( 11 . 1 . 1 ) 


where and are the previously derived tangent stiffness matrices for 

the superconducting and thermal conduction problems respectively and 
is the pre^dously derived loading vector for the thermal conduction problem 
with convection boundary conditions. The matrix is df^^fdv'^ and 

the zero entries on the left euid right hand sides of the equality sign appear 
because /dv^^ = 0 and = 0. Because the resistivity of each 

element is zero for a superconductor, only one nonzero term appears in 
q^. This term has a magnitude of 7 l ^md appears at the degree of freedom 
for T associated with the radial distance r^. 

Solving the set of thermal equations produces the simple result that 
T is a constant over the whole conductor, the value of T at each node being 
equal to 7^ + Cn+i^- This simplifies the solution procedure considerably 
because the thermal equations = q^, need not be solved by assem- 

bling and inverting K^. Instead of solving (11.1.1), which also requires the 
assembly and storage of the corrective Newton-Raphson technique 

is used to directly iterate to a solution for This is accomplished by 

solving the superconductor problem with the current load I held constant. 
The value for T at each node is the value of T at the n + 1 step, 7^ + Cn+i^ » 
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where Cn+i equeils Cn + and the quantity /„ is the input step size. This is 
the value of T that is also used to determine a and since they are both 
functions of T and their values cire required to obtain the correct EM tangent 
stiffness matrix 

Essentially this is the same incremental/iterative solution method 
used in the two previous chapters to generate the exact solutions for 

and The loading parameters and axe held constant and 

the system is allowed to iterate onto the equilibrium path. As mentioned in 
those chapters, this technique can fail or become computationally expensive 
when a Izirge niimber of corrective iterations is required. Early in the testing 
of the STEPlD element, it was observed that a reasonable solution for 
and could be obtained in this manner if the reference state is set so that 
all unconstrained values of are equal to iV’ool values of are 

set equal to zero. The niimber of iterations needed for convergence with 
this reference configuration was usually less than ten. This configuration is 
also identical to the initial reference state that is used to start the solution 
procedure for a superconductor when the current loading I is varied. This 
means that the subroutine that generates the initial reference state for the 
varying current problem can also be used to generate a reference state for 
the varying T problem thereby reducing the logic and memory requirements 
of a code that solves both problems. 

The computational cost of using the above solution method is that 
more iterations are required at each step and the reference configuration must 
be recomputed at each step. There is also no guarantee that the solution 
method will convege but numerical experiments strongly suggest that it will. 
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The corrective Newton-Raphson technique is therefore used to solve 
the problem of a superconductor with a varying thermal boundary load for 
the following reasons: 

1. Assembly and storage of is unnecessary. 

2. The number of iterations required per step is relatively modest. 

3. Convergence, as determined by numerical experiments, always occured. 

4. The reference state is identical to the initial configuration of the I load- 
ing problem, allowing the use of one subroutine to set either configura- 
tion. 

11.2 DETERMINATION OF THE CORRECT EQUILIBRIUM PATH. 

To check whether the conductor is in the normal or superconducting 
state, one determines the critical temperature "Tc and the critical magnetic 
field Be- Then and T; are compared to the largest magnitude of B and the 
highest temperature T (typically Too) field within the conductor. If either T 
or Be is exceeded, a superconductor changes quantum states and becomes a 
normal conductor. This change of state occurs because at T equal to T or B 
equal to B. a bifurcation point for the equilibrium paths of a superconductor 
and a normal conductor exists. 

The existence of this bifurcation point also means that if the con- 
ductor is originally in the normal state and Too falls below T o,ni the largest 
value of B within the conductor is less than B^, the material becomes su- 
perconducting. It is therefore important to know the values of 7^ and Be so 
that the position of the bifurcation point along the equilibrium path can be 
determined. The critical temperature 7^ is a material constant, and is de- 
termined either by experimentation or by referencing previous experimental 
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data. The critical B field Be is determined by using Equation (4.4.1) which 
is a function of the temperature T. 

An alternative method for finding the correct conductor state is 
to compute the Helmholtz free energy for the superconducting and normal 
states of a conductor for the same thermal and current loading conditions. 
After finding the free energies of both systems, the state of the system can 
be determined by choosing the system that has the lower free energy. This 
approach is computationally inefficient because it requires solving for the de- 
grees of freedom v of both states at every incremental step. It also requires 
knowledge of the heat capacity of the material for the superconducting and 
normal states of the material, and finding these values can be a task of con- 
siderable difficulty because of the dearth of experimental data. The first 
approach is therefore chosen here. ^ 

Following the first approach, the B fields and T distribution as de- 
termined by axe checked at the end of each incremental step to see 

if they axe sufficiently small so that a superconducting state is possible. If 
that is the case and the system was originally in the normal state at step 
n, then v„+i is solved for again at step n + 1 using the superconducting 
finite element STEPlD. If that is not the case and the condutor was in the 
normal state at step n, then the solution at step n -f- 1 is accepted and the 
program proceeds to the next step keeping the normal conducting state ele- 
ment LET ID. If the B fields and the T distribution axe small enough at step 
n -f 1 and the system was originally in the superconducting state at step n, 
the solution vector v„+i is accepted and the solution process moves to the 
next step using the STEPlD finite element. Finally, if either of the two fields 
axe too laxge for a superconducting state to exist, then is recomputed 
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using the LETID finite element snd the solution method moves to the next 
step using the LETlD finite element. 

For the one-dimensional problem, this methodology appears to be 
optimal and is the one used for the present numerical experiments. It has 
the key advantage that knowledge of the heat capacity of the material is 
unnecessary. Furthermore, the solution vector v only needs to be computed 
twice when the system changes state. If the second path determination 
method that involves computing the Helmholtz free energies of the normal 
and superconducting states is used, the solution vector v must be computed 
twhce at every step. 

The first approach naturally delineates the tests for the proper equi- 
librium state of the system into four separate cases where a change of state 
may occur. These cases are: 

1. System originally in the superconducting state, thermzil load increasing. 

2. Syst em originally in the superconducting state, current load increasing. 

3. System originally in the normal state, thermal load decreasing. 

4. System originally in the normal state, current load decreasing. 

For cases where the system is originally in the superconducting state 
and the current or thermal loading is decreasing, the system remains in 
the superconducting state because the solution is moving away from the 
bifurcation point and no problems involving equilibrium path changing are 
posed. Similarly, when the system is originally in the normal state and the 
current or thermal loading is increasing, the solution remains in the normal 
state because it is moving along the normal state equilibrium path away 
from the bifurcation point. Again, this poses no problems to the solution 
procedure. 
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The comparison of % to Too, the first path changing criterion, is rela- 
tively straightforward and only involves the computation of 7^ at each step. 
The second criterion can cause problems. This criterion states that B for 
every element of the conductor must be below Be for a superconducting state 
to exist. Comparing B over each element to Be can become computationally 
expensive as the number of elements used to model a problem increases. This 
computational expense can be reduced by finding a priori where B attains 
its largest %'alue within the conductor. For the one-dimensional axisymmetric 
infinite conductor, B attains its largest value at rg. A cursory examination 
of Equation (5.1.8) verifies that the preceding statement is correct. Equation 
(5.1.8) can also be used to determine that 

^ ( 11 . 2 . 1 ) 

where /x<> replaces for example conductors of this work, as discussed in pre- 
vious chapters, and j • ficdT is equal to I. By using this analytical solution 
for Be(rc). a simple means exists for determining if a superconducting state 
is possible. 

The only situation not discussed so far is when the incremental solu- 
tion falls directly upon the bifurcation point. As explained previously, at this 
point the physical solution cannot be modeled by one-dimensional elements 
and the LETID element is used to model EM quantities although the solu- 
tion generated does not represent the correct physical state of the system. 
This method is used so that the solution method can proceed to the next 
step without failing. The LETID element does not fail at this point because 
it is based upon a potential energy formulation. The STEPID element fails 
because it is based upon the difference of the Helmholtz free energies of the 
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superconducting and normal states. At the bifurcation point, this difference 
is zero and the tangent stiffness matrix becomes singular. 

The physical significance of what is occuring is that both the normal 
and superconducting states possess the same energy. In the actual physics, 
a variation of ^ occurs in the z direction, and the system chooses the eigen- 
state that possesses the lowest possible energy and entropy. To extend the 
current superconducting model to this problem, an adaptive mesh appears 
to be necessary to determine the boundary between parts of the conductor 
that are normal and superconducting. The adaptive mesh is also required to 
make well-conditioned enough that reasonable values of EM quantities _ 

can be generated by standard nonlinear solution techniques. Unfortrmately, 
time limitations on the thesis research precluded the development of an adap- 
tive two-dimensional mesh and the examination of the physics of this most 
interesting and challenging problem was foregone. 

The STEPID finite element used for the thesis research should pos- 
sess a rank deficiency of one at the bifurcation point. In an effort to gain 
a better understanding of what was occuring at the bifiurcation point with 
the finite element model, the model was forced to converge upon this point 
by setting the thermal loading to the critical temperature % and setting 
the current loading to a value that would generate the critical field Be at 
the outer conductor boundary The corrective Newton-Raphson solution 
method was then used to iterate to the bifurcation point. The finite element 
model actually converged and returned a quantum state that carried no cur- 
rent and an applied external field of Be at r^. Even though should be 

singular at this point, a fact that precludes convergence, it is believed that 
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the STEP ID model converged for two reasons. First, the CNR iterative pro- 
cedure is stopped when the 2-norm of r is smaller than the input tolerance 
r. Second, the sccilings and the factorization of introduce numerical 

round off errors that perturb the generated solution off of the bifurcation 
point just enough to render nonsingular. 

The STEP ID model returned the restilt of an applied external field 
and no current in the conductor because it does not enforce the current 
conservation constraint I — fp dTAc • j = 0. The addition of this constraint 
automatically allows an EM model to distinguish between cases where the 
field at Tc is generated by a current / or by an externally applied B field. 
The eigenvalue analysis of for earlier versions of Ginzburg-Landau and 
London superconducting finite elements that contained the current conserva- 
tion constraint showed that the current conservation constraint is redundant 
when no external fields are applied to the system or when does not lie 
upon the bifurcation point. These two cases aire not considered in this work 
and axe not presented here. 

To summarize, the basic path determination process is as follows: 

A. Solve system of equations for 

B. Find Too by using Zx, = To -h Cj+i'Ti- 

C. Find Be at T» by using Equation (4.4.1). 

D. Find I by using I = Io + C^^+i-^L- 

E. Find B$ at Xc by using Equation (11.2.1). 

F. If B$(rc) > Be, go to H.; if not, go to G. 

G. Toe >Tc, go to H.; if not, go to I. 

H. If the current EM element type is STEPID, change it to LETID and go 
back to A.; if it is not, go to J. 
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I, K the current EM element type is LETID, change it to STEPID, reset 
the reference state of and go back to A.; if it is not, go to J. 

J. Accept solution vector v„^-i for the current step. 

Additional logic statements are included in the actual coding to help 
prevent the STEPID element from exceeding the bifurcation point. Beyond 
the bifurcation point, the element will “zero in” upon the same type of solu- 
tion that occurs when the element is forced to converge upon the bifurcation 
point. A quantum state is generated where there is no current in the con- 
ductor and a boundary B$ field loading of magmtude exists. 

This state usueJly requires more iterations to converge upon a solution and a - 
larger solution tolerance r th 2 in physical states that lie below the bifurcation 
point on the equilibrium path. Non-physical solutions for the superconduc- 
tor that lie beyond the bifurcation point can also cause the CNR procediire 
to fail if T is too small or if the maximum number of allowed iterations for 
the solution procedure is exceeded. 

To prevent the STEPID element from moving past the bifurcation 
point and encoimtering these problems when the current load I is being 
incremented, steps C through E of the path determination procedure are 
performed prior to each corrector iteration. The current iteration value of 
C^^+i is used for step D because the actual step size along the eqtiilibrium 
path may change with each iteration. If the step size changes, then 
changes for each iteration and so does the value of /. If Be is exceeded at Tc 
for any iteration, the program changes the element type to LETlD, and the 
solution procedure restarts at step n and attempts to increment the loading 
to step n + 1. 


Ililllilillll I 
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For the case where the thermal loading is varied, steps B through 
E are performed prior to each predictor step. For this case where the tem- 
perature T is being incremented, the current 1 is held steady and causes no 
problems because is held constant and is known before the predictor and 
corrector steps are taken. T is also known a priori as being 7^ -I- (Cn "b 
at step n -F 1 as discussed in Section 11.1. Because Cn+i known before 
the solution procedure begins, steps B through E of the path determination 
process can be used to determine if the solution vector will move past the 
bifurcation point before the solution process begins. Inserting this test before 
the predictor step keeps the program from performing an unnecessary- solu- 
tion step. After steps B through E of the path determination process axe 
performed, Bc(T^-i-i) and Tc are compared to B${rc) and Tn+i respectively. 
K either of the two latter quantities exceed their critical values the EM el- 
ement type is changed to LETID and the solution procedure is allowed to 
continue. If the critical values axe not exceeded, the solution procedure is 
allowed to continue imaiFected. 

11.3 NUMERICAL EXAMPLES. 

The LINTID, LETID and STEPID finite elements derived in previ- 
ous chapters have been applied to the solution of the test problems described 
later in this section. The CUPLEID and STEPID elements were used to 
determine EM quantities and the LINTID element was used to determine 
the thermal distribution for the normal state of the superconductor. The 
temperature of each node of the superconductor is calculated as described in 
Section 11.1 of this chapter. The description of the nodal degrees of freedom 
and the calculation of the material properties of each element may be found 
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in their respective chapters. The application o£ boundary conditions, scal- 
ing techniques, mesh generation techniques and the assembly and solution 
techniques are also described in the respective chapters. For the graphical re- 
sults generated in this chapter, B$ was calculated by using the integral form 
of Maxwell’s inhomogeneous equation for magnetic fields, Equation (7.1.5). 
This equation requires knowledge of to determine Be. For the super- 
conducting phase, is calculated by using Equation (3.2.16) and for the 
normal phase is calculated by using the elemental values returned in the 
solution vector v. The integral of Equation (7.1.5) is evaluated by 2 point 
Gaussian integration. 

11.3.1 PT?ORT.KM 1: VARYING T LOAD. 

For this problem and the next, the test material is high purity alu- 
minum. Reference and initial states of the system are set as described in 
previous chapters. The geometry is that of a one-dimensional axisymmetric 
wire as shown in Figure 2.1. The wire radius Tc is 1.15 x 10“^ meters and 
transported a total current I in the positive z direction. The initial current 
lo is 1 ampere and the loading current II is 0 amperes. The initial free 
stream temperature is .5625° Kelvin and the loading temperature is 1 
Kelvin. Because the free-space magnetic field element has been validated 
previouslv. all elements are within the wire. The mesh for the superconduct- 
ing phase has 98 elements in the boundary layer and 2 elements in the bulk 
layer while the mesh for the normal phase of the conductor has 50 elements 
in the coarse mesh and 50 elements in the fine mesh. As in Chapters IX 
and X the depth of the fine mesh was .25 Vc for the normal conductor. The 
depth of the boundary layer mesh for the superconducting phase varies with 
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temperature as discussed in Chapter VIII. The step size /„ is chosen as .0125 
and 80 steps eure used to increment from 0.0 to 1.0. A solution tolerance 
r of 4 X 10“^^ is used for the superconducting state cind 9 x 10~^ is used for 
the normal state. 

The solution procedure required 41 steps in the superconducting 
phcise averaging 4.61 iterations per step. The estimated condition number 
for these steps ranges between 181 and 834. The solution procedure then 
required 39 steps in the normal phase with an average of 2 iterations per 
step. The estimated condition number varies between 29861 and 204664. 

Data output files for all of the figures to be shown for all examples 
in this chapter axe saved every tenth step and all graphical representations 
of this data are labeled with the appropriate values of or when 

it was possible. Graphical representations of each data set axe generated 
by using the PLOT2D utility to produce a raster file and then using the 
raster files to create a PostScript Izmguage file. This is mentioned because 
the graphical representations of data sets axe subject to the limits of the 
PLOT2D utility. The data sets for each variable were then loaded into a 
single file, and PostScript language commands were used to generate axes 
and data set labels and legends. This is mentioned for researchers who 
wish to duplicate the graphical results because the PLOT2D utility does not 
possess the ability to add the desired labels and font types or graph 10 sets 
of data on a single graph. 

Restdts for the temperature distribution within the wire axe shown 
in Figure 11.1 and match the expected physical behavior. The results for 
IV’l^/lV’ool^ in the region 1.023 x 10“'* ^ r < are presented in Figure 
11.2. The value of the normalized value of 10| is not shown over the whole 
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conductor because all of the physics of interest occured within the boundary 
layer. The value for iV’lVl’/’ool^ for r in the region 0 < r < 1.023 x 10~'‘ is 
a constant equal to 1.0. The expected physical behavior for the normalized 
value of 1^1 in the boundary layer was that as the temperature increased, 
the boundary layer depth would increase and |V’| /IV’ooP would vary over 
a wider range of r. This physical behavior is accurately captured by the 
STEPID element and is shown in Figure 11.2. Figures 11.3 and 11.4 show the 
value of the current density in the normal and superconducting phases 
respectively. For the superconducting phase, only the boimdary layer values 
are shown with all other values of being equal to zero. The value of jr ^ 
at each node for the superconducting phase is calculated by use of Equation 
(3.2.16). Because the current I was steady, the magmtude of the current 
density should decrease as the temperature increases for the superconducting 
phase. The boundary layer depth shoxild also increase. Again both physical 
characteristics are accurately depicted by the STEPID model. 

For the normal phase of the conductor, is depicted as a step 
function in Figure 11.4. The step function represention is necessary because 
the current density is approximated by a step function over an element by 
the LETID finite element. The results for the steps where equals 0.7 
and 0.9 were omitted for clarity. By referring back to Figure 11.1 some 
determinations can be made about the behavior of ^ for the normal phase. 
It can be seen that the temperature is higher at the center of the conductor 
than at Vc- The resistivity should also be higher at the center and ji ^ should 


be smaller there. 
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As the temperature of the wire increases, the amount of thermal 
energy produced by the current I through a wire should remain almost con- 
stant. The rest of the thermal energy in this system comes from the free 
stream temperature boimdeury conditions. The significance of this is that 
the temperature distribution in the wrire should become more homogeneous 
and the magnitude of the heat absorbed by the system from the free stream 
boundary conditions should eventually become greater than the magmtude 
of the heat produced by resistance to the current I. The resistivity will also 
be determined more by the boundary conditions than by the heat generated 
by the steady current I. The temperature should become more homogeneous 
throughout the wire and the resistivities and current densities should follow 
suit. This expected behavior is accurately modeled by the fimte element 
approximation as can be observed in Figure 11.4. 

The only behavior that at first appears to be non-physical is the 
jump in the magnitude of the current density as the conductor changes from 
the normal to the superconducting phase. This is easily explained because 
is a function of the resistivity a; and the E field for the normal state while 
it is a function of and \tp\ for the superconducting state. The easiest 
way to verify that the current density predicted by the finite element method 
is exact is the use of the integral form of Maxwell’s inhomogeneous magnetic 
field equation (Equation (7.1.5)) to evaluate Bg at rg. This equation requires 
that no matter what the jz distribution may be, that for wires carrying the 
same current 7, the value of Bg at Tc will always be the same. Because the 
current I is held steady for this example, Bg should eJways be the same at Tc 
independent of the quantum state of the conductor. Figure 11.5 shows this 
expected behavior accurately and also demonstrates why Equation (7.1.5) 
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has been used to compute the Bg field for the results presentation of this 
chapter. Equation (7.1.5) aUows the value of Bg to be computed at each 
node while the equation used in previous chapters, Equation (5.3.2), only 
allows the computation of the mean value of Bg over each element. By 
using Equation (7.1.5), the accuracy of computed values of can easily be 

verified by comparing values of Be at Tc- 

Figures 11.5 and 11.6 show the distribution of the Bg field within the 
conductor as the temperature was increased. Figure 11.5 shows only values 
of Bg that lie within the boundary layer while Figure 11.6 shows values for 
the normal state of the conductor for r between 0 and Tc. In Figure 11.5, it 
can be seen that the Bg field penetrates more deeply into the conductor as 
the temperature of the conductor increases. This is the desired and expected 
physical behavior. In Figure 11.6, it can be seen that the small increase in 
the temperature for equal to .6 to 1.0 produces no significant changes 
in the Bg field. The changes are so small that the PLOT2D utility connot 
discern changes in Bg as the temperature increased although there is a small 
change in the Bg field that follows changes in The maximum nodal 

change of Bg as is varied from .6 to 1.0 was ~ 2 x 10““ percent. This 
percent difference of Bg occurred between the states where was equal to 
.6 and 1.0. 

11.3.2 PROBLEM 2: VARYING I LOAD. 

As mentioned earher, the test material is high purity aluminum. 
Reference and initial states are set as described in previous chapters. The 
geometry is the same as that of Problem 1 of this chapter with being 
equal to 1.15 x lO"'^ meters, Ncoarse = Nfine = 50 elements for the normal 
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state and Nbuik aJid Nbound being equal to 2 and 98 elements respectively 
for the superconducting state. The initial current !„ is 1 ampere and the 
loading current II is 2.1 amperes. The initial free stream temperature % is 
1° Kelvin and the loading temperatme 7 l is 0® Kelvin. Again no elements 
were generated external to the conductor/free space boimdaxy located at Tc. 
Meshes for both states axe generated as described in Chapters VIII and IX. 
The step size l„ is chosen to be .01 and 80 steps were used to increment 
from 0.0 to .73. The solution tolerances r of 4 x 10"^’’ and 9 x 10“^ are used 
for the superconducting and normal states respectively. 

The solution procedure required 40 steps in the superconducting 
phtise averaging 3.65 iterations per step. The estimated condition number 
for these steps ranged between 231 and 10935. The solution procedure then 
required 40 steps in the normal phase averaging 2.05 iterations per step. 
The estimated K condition number for these steps varied between 28146 
to 181319. Data output files were saved every tenth step as stated in the 
previous section. 

Results for the temperature distribution in the whole wire axe shown 
in Figure 11.7. Because the free stream temperature is held constant at 1° 
Kelvin, the major source of heat energy comes from resistance of the current 
fiow I through the wire rather than from boimdaxy loading. This caused the 
temperature differential between the center of the wire and r equEil to Xg to 
be greater than in Problem 1 of the previous section. This expected physical 
behavior is shown in Figure 11.7. 

Figure 11.8 depicts the behavior of IV’l^/lV’ooP Ibr the region where 
r varies between 1.0597 x 10~^ eind 1.1887 x 10~^ meters. For r betvreen 0 
and 1.0597 x 10“^ meters, |i/’|^/lV’ooP is unity. Graphical results of the data 
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obtained for equal to .16 and .30 are omitted for clarity. Physically 
it is expected that as the current I is increased, the system will move to a 
higher energy state. As the energv^ of the system increases, the botmdary 
layer should widen independently of whether the energy sotirce is thermal 
or electromagnetic in nature. This behavior is accurately reflected in Figure 
11.8, but comparison of this plot with Figure 11.2 shows that for the current I 
loading case, it appezirs that there was more of a shifting of the distribution of 
the Cooper pairs towards the center of the conductor rather than a reordering 
of the distribution. An explanation of the physics of these two cases can be 
made by invoking the London model of superconductivity. 

2 

With the London model, the number density of Cooper pairs |V’| 
is only a function of the temperature T and is equal, to |r/»ooP- Using this 
information and refering to Equation (4.4.4), it can be seen that as the 
temperature increases, the total number of the Cooper pairs will decrease. 
This is the general behavior of the Ginzburg-Landau superconductor as well. 
As the temperature increases in Figure 11.2, the number of Cooper pairs 
in the current stream must also remain constant because I is constant. To 
maintain the same number of Cooper pairs within the current stream as ”7” 
increases, the system must reorder itself and impart a kinetic momentum to 
pairs that lie deeper within the boundary layer. 

For the case of an increasing current I and steady temperature T , 
the total number and distribution of the Cooper pairs must remain approx- 
imately constant. As the current I is increased, the number of the Cooper 
pairs remains essentially constant, but the number of Cooper pairs with a 
kinetic momentum increases. This increase in the energy of the system de- 
stroys some of the Cooper pairs. The annihilation of Cooper pairs occurs at 
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the most energetically favorable position, within the boundary layer. This 
also serves the dual purpose of widening the boundary layer, upon which 
an increased number of Cooper pairs with kinetic momentum is allowed to 
move more easily. This expected behavior explains why there is more of 
a “shifting” of the distribution of the Cooper pairs in Figure 11.8 than a 
reordering of the distribution. The boundary layer is widening in response 
to the increasing number of pairs with a kinetic momentum. In Problem 1, 
the amount of Cooper pairs annihilated by the increasing thermal energy is 
much greater than those destroyed by the increasing EM energy of Problem 
2. This relative change in the number of Cooper pairs as an incremental step 
is taken explains the natiire of the difference of the two plots. 

Figures 11.9 and 11.10 show the change in as was incre- 
mented. Figure 11.9 shows the superconducting state emd Figure 11.10 shows 
the normal state. Figure 11.9 only shows values in the boundary layer. Out- 
side of this layer, the plotted values of for the superconducting state 
vanish. For this figure, different line types and a legend are used so that 
the plots for equal to .30 and .32 are more easily distinguishable. The 
graphical representation of in Figure 11.9 matches the expected physi- 
cal behavior. As I is increased, the boundary layer should spread and the 
magnitude of should also increase. The STEP ID finite element also 
captured this expected behavior well. As in Problem 1, the current density 
should be higher at Tc than the center of the conductor for the normal state. 
There should also be an increase in the magnitude of 8is the current I 
is increased. The LET ID element performed as expected and modeled this 
expected physical behavior. 
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Figures 11.11 and 11.12 show the Finite Element values of Be plotted 
over the boundary layer and the whole conductor respectively. Figure 11.11 
omits the labels for equal to .50 and .66 for clarity. Similarly, Figure 
11.12 omits the labels for equal to .16, .23, .30, .50 and .66. The ex- 
pected physical behavior for the superconducting state is that as the current 
increases, the magmtude of Be will increase and penetrate more deeply into 
the boimdary layer. For the normal state, the magmtude of Be should keep 
on increasing but it should also be an almost linear function of the distance 
r from the center of the conductor. Both of these physical behaviors are 
again modeled well by the finite element computed solutions and illustrate 
the ability of the four-potential based finite elements to model EM fields. 
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Figure 11.8: IV’lViV’ooP vs- r. Values for |t/>|VlV’ooP 1 plotted. Graphs for = .16 and .30 omitted for clarity. 
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Figure 11.11; Bg vs. r. Nonzero values for the superconducting state plotted. Values for the normal state also plotted 
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11.4 SUMMARY. 

This chapter “glues” together all of the previously derived finite el- 
ements into a comprehensive program that can determine the correct equi- 
librium state of a thermally and quantum mechanically coupled EM system. 
The primary emphasis is the discussion of the results for two different cou- 
pled problems but two other topics are also discussed: the solution of the 
superconducting problem where I is constant and Too is varied, emd the de- 
termination of whether an EM system is in the superconducting or normal 
state. 

The constant I, varying Too problem mentioned above is solved 
rather easily 8is is the determination of the correct quantum state. The 
only real problem and failing of the final model that is developed here is its 
inability to accurately model a system within a conductor that has mixed 
normal and superconducting states near the transition point. Fortunately, 
the four-potential method is readily extensible to the solution of this problem 
although it is not addressed in this thesis primarily because the development 
of an adaptive two-dimensional mesh to deal with conditioning problems of 
the tangent stiffness matrix would have required a considerable investment 
of time and effort. 




CHAPTER XII 


CONCLUSIONS. 

12.1 SUMMARY OF WORK. 

As mentioned in Chapter 1, the primary purpose of this work is 
to develop a finite element model for types I and II superconductivity that 
can accurately predict EM quantities. This model is to include thermal ef- 
fects and to have the ability to change between superconducting and normal 
phases when necessEiry. Originally, this model was to be based upon the four- 
potential variational principle to reduce the number of degrees of freedom 
per element node. However, it weis discovered during the course of research 
that the four-potential variational principle offered more advantages for the 
analysis of EM problems than just the simple reduction of element nodal de- 
grees of freedom. More important is the ability of the four-potential method 
to model any EM problem that heis been posed here through the adjimction 
of constraints by a Lagrangian multiplier. An equally important advantage 
of the four-potential method is that B and D discontinuities at material 
interfaces cire enforced automaticadly and require no special attention from 
the user. The current predicting elements presented in this work required a 
special boundary treatment solely because j is used as an independent vari- 
able instead of $. This choice weis made originally to simplify the cmrent 
conservation constraint and was not changed. The simpler $ formulation 
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given by Equation (3.1.5) reduces the number of degrees of freedom required 
by a one-dimensional current predicting finite element by two. 

To produce the desired four-potential based EM finite element men- 
tioned in the first paragraph, a functional that used A and $ is first de- 
veloped for any arbitrary material. This functional is then augmented by 
the Lorentz gauge constraint to ensure that A is unique. The augmented 
functional is then applied to one and two-dimensional geometries and the 
natural boundary conditions of the two geometries are determined. At this 
point, it was determined that a further extension of the new functional was 
necessary to model the more general case of an unknown current density j. 
This extension is necessary because the arbitrary nature of geometries for 
EM problems does not always permit a priori knowledge of the distribution 
of a current within a conducting medium. It was also realized that tempera- 
ture differentials within a conductor make the resistivity within a conductor 
inhomogeneous. The varying resistivities also preclude an a priori knowledge 
of j within a conductor. 

To model the thermal effects that are eventually added to the EM 
model of a normal state conductor it was therefore necessary to extend the 
previouslv derived four-potpntial functionals to include cases where j is un- 
known. This is accomplished by augmenting a gauged form of the four- 
potential functional by an additional constraint, the current conservation 
constraint. The function^ is also modified by making j a primary variable 
instead of # through the use of the constitutive relation between $ and j. 
This substitution requires the additional augmentation of the functional by 
a boimdary continuity constraint. The additional constraint is necessary 
because the previous substitution for $ inhibits a necessary integration by 


236 


parts that ensures the continuity of the E field across material interfaces. 
The final functional is a four-potential based functional for determining EM 
fields in hneax conducting materials. 

The next phase extends the four-p>otential variationed principle to 
cover type I and II superconducting materials. The Ginsburg-Landau equa- 
tions pro\'ide the necessary basis for this extension. The variational func- 
tioned used to derive these equations contains the magnetic vector potential 
A. This functional also contains terms that represent a Landau expansion of 
the Helmholtz free energy of the quantum wave order peirameter ^ around 
the critical temperature 7^. To adapt this functional to the four-potential 
method the electric field energy Ue is added and the gauge constraint ad- 
joined. The gauge constraint used here is the London gauge which is identical 
to the Lorentz gauge for magnetostatic problems. Because all of the super- 
conductivity cases that are considered here are free of electrostatic charge, 
the electric field energy Ue is zero and this term is not included in the 
functionals of this work. After the augmentation of the Ginzburg-Landau 
variational functional by the London gauge conscraint is complete, the two 
material parameters a and /3 of the Landau expansion are determined as 
functions of the effective penetration depth A^// and the critical magnetic 
field Be of a superconducting material. 

Conventional thermal field variational functionals are then used to 
describe the thermal energy of the EM systems imder consideration. The 
temperature dependence of material parameters is also developed for con- 
ductors in both the superconducting and normal states. 

This modeling work completes the necessary background for the de- 
velopment of EM finite elements that are thermally and quanttim phase 
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coupled. Elements and solution procedures are then developed. The most 

significant features of the finite elements are: 

1. The normal element has the ability to predict current densities with a 
high degree of accuracy. 

2. The superconducting element has the ability to show the current density 
distribution in much greater detail than ever before. The significance of 
this feature is that if the Ginzburg-Landau model of superconductivity 
is correct, there is a greater imderstanding of the physics that occur 
within a superconductor. 

3. A nonlinear superconducting finite element that does not require path 
follovdng procedures to determine equilibrium states if the correct ref-, 
erence state and mesh are chosen. 

4. A superconducting finite element that also is rapidly convergent upon 
a solution, well conditioned and, as far as the author has been able to 
determine, generally convergent upon the equilibrium solution provided 
the correct reference state and mesh have been selected. 

5. The combined use of the thermal, normal, and superconducting el- 
ements provides for a comprehensive program that can analyze any 
physical eqmlibrium state of a conductor except for the mixed nor- 
mal/superconducting' state. Appropriate modifictions to allow for the 
modeling of this state are also suggested. 

6. Finite element models that can model any EM media provided that the 
thermal and EM properties of the medium are known. 

7. Finite element models that are modular and employ standard linear and 
nonlinear assembly, scaling, and solution techmques. 
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8. Finite element models that require no special boundary treatment for 
adjacent elements that possess differing EM or thermal properties. 

9. EM finite elements that can predict electric and magnetic fields with a 
high degree of accuracy. 

10. EM Finite elements that require fewer degrees of freedom for the analysis 
of two and three-dimensional field problems than the conventional field 
based finite elements currently in use. 

12.2 DIRECTIONS FOR FURTHER RESEARCH. 

The focus of this work is upon the analysis of magnetostatic EM field 
problems. These cover a significant but small part of the range of EM field 
problems that axe of interest to scientists and engineers. The ready exten- 
sion of the four-potential variational principle to a wide range of EM field 
problems provides a powerful tool for the solution of difficiolt EM problems. 

An unsolved problem of most interest to the author is the one where 
normal eind superconducting state coexist near the transition state. A re- 
alistic treatment of this problem requires two- and three-dimensional space 
discretization and consequently follows outside the scope of this work. To 
extend the present work to that problem, a multidimensional adaptive mesh 
appears to be necessary to determine the interface between normal and su- 
perconducting portions of a conductor as well as to improve the conditioning 
of the tangent stiffness matrix. 

Another direction that scientific research can take from the results of 
this work is experimental verification of these results. This verification would 
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add considerably more evidence for the validity of the Ginzburg- Landau the- 
ory of superconductivity as well as a greater understanding of supercon- 
ductivity in general. Even without this verification, the results herein that 
compare the current density to a low viscosity fluid stream might be applied 
to the aneJysis and eventual development of a model of high temperature 
superconductivity. 

Other problems of active interest to the engineering community are 
dynamical in nature. These problems include the analysis of time- dependent 
EM waves moving through fixed (static) EM media, as well as EM media 
coupled with rapid mechanical motions. These problems are highly comples, 
but the general applicability of the four-potential method to EM problems 
in general appears to be well suited for the numerical treatment of these 

problems. ■ :n. -q r 

Finally, the thermal functionals that are used to analyze the tem- 
perature distribution within the conductor are adequate for the relatively 
minor loadings and changes of loadings that are presented here. A direction 
of fiirther research that the author has already undertaken is the develop- 
ment of thermal finite elements that are nonlinear in nature. These elements 
allow the thermal conductivity k and the electrical resistivity w to be func 
tions of the temperature T rather than the spatial coordinates. It is hoped 
that these elements will permit a more accurate analysis of the temperature 
distribution within a normal conductor and allow larger steps to be taken 

with solution path following techniques. 

In conclusion, the ELI finite elements that are presented here have 
performed well and confirmed the ability of the four-potential variational 
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principle to solve a range of problems. The author believes that this method- 
olgy is relatively simple to use and exhibits key advantages over current field 
based formulations. Potential based formulations and variational principles 
show promise for the treatment of unsolved EM problems and should both 
be given due consideration over Seld-b€Lsed formulations. 
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